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A  wide  variety  of  problems  in  different^!  and  integral  equations  require  ab- 
plication  and  inversion  of  linear  operators.  For  large-scale  physical  problems,  the 
size  of  the  problem  n  is  such  that  the  work  esmended  by  general  algorithms-0(n*) 
for  application  of  a  transformation  and  O(n^)  for  application  of  its  inverse-is  of¬ 
ten  prohibitive.  Several  methods  have  been  devised  in  recent  years  to  reduce 
these  complexities  by  exploiting  the  structure  of  particular  problems. 

This  thesis  puts  earlier  methods  into  a  unified  framework  by  developing  new 
wavelet-like  bases  for  /2^[0, 1],  and  other" spaces,  which  lead  to  efficient  algorithms 
for  operator  application  and  inversion.  It  is  based  on  the  observation  that  in 
many  cases  of  interest,  while  the  matrices  involved  are  dense,  their  elements  vary 
smoothly  as  a  function  of  their  indices,  except  along  a  collection  of  bands  of  fix^ 
width.  Algorithms  are  presented  for  transforming  such  locally  smooth  map*ftes 
into  matrices  which  are  sparse  (to  a  finite  but  arbitrarily  high  accurapy^  The 
resulting  sparse  representations  support  fast  algorithms  for  matrix^pplication 
and  inversion,  requiring  variously  0(n),  O(niogn),  and  O(nlog^)  operations. 

Programs  have  been  developed  applying  these  techniques  to  the  solution  of 
second-kind  integral  equations  with  non-oscillatory  kernels.  Numerical  results 
are  presented  to  demonstrate  the  effectiveness  of  the  approach.  '  r  , 
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Chapter  1 
Introduction 


1.1  Background 

Integral  equations  are^  a  well-known  mathematical  tool  for  formulating  physical 
problems.  Historically  they  have  achieved  great  popularity  among  mathemati- 
\  cians  and  physicists  in  formalizing  boundary-value  problems  of  gravitation,  elec¬ 
trostatics,  fluid  dvnamics,  and  scattering.  Another  application,  of  considerable 
current  interest,  is  semiconductor  modelling.  While  a  number  of  physical  prob¬ 
lems  give  rise  to  integral  equation  formulations  directly,  a  wide  variety  of  partial 
differential  equations  may  also  be  expressed  as  integral  equations. 

Integral  equation  formulations  have  several  advantages  (good  conditioning, 
dimensionality  reduction,  and  the  ability  to  treat  arbitrary  regions)  but  have 
had  one  overriding  disadvantage:  the  high  cost  of  working  with  the  associated 
dense  matrices.  For  a  problem  requiring  an  n-point  discretization,  the  inverse 
of  a  dense  n  x  n-matrix  must  be  applied  to  a  vector.  Even  to  apply  the  matrix 
itself  to  a  vector  requires  order  O(n^)  operations;  application  of  its  inverse  by 
a  direct  (non-iterative)  method  requires  order  O(n^)  operations.  If  an  iterative 
method  is  employed,  the  number  of  iterations  depends  on  the  condition  number 
of  the  problem  and  each  iteration  requires  application  of  the  n  x  n  matrix.  For 
large-scale  problems,  the  resulting  costs  are  often  prohibitive. 

In  recent  years  a  number  of  algorithms  ([2],  [8],  [9],  [15])  have  been  developed 
for  the  fast  application  of  linear  operators  naturally  expressible  as  dense  matrices, 
the  best  known  of  which  are  the  particle  simulation  algorithms  developed  by 
L.  Greengard  and  V.  Rokhlin  [8].  Each  algorithm  of  this  class  exploits  the  special 
structure  of  a  particular  problem  by  combining 

1.  interpolation  of  the  function  which  defines  the  matrix  elements,  with 

2.  a  divide-and-conquer  strategy. 

For  example,  the  particle  simulation  algorithms  are  based  on  the  observation  that 
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the  electrostatic  or  gravitational  field  due  to  a  set  of  particles  varies  smoothly  in 
regions  separated  from  the  location  of  the  particles  (see  Fig.  1.1).  Such  a  field 
can  therefore  be  represented  to  high  precision  by  an  expansion  containing  only  a 
few  terms.  This  observation  is  augmented  by  the  construction  of  a  divide-and- 
conquer  scheme  in  which  virtually  all  particle  interactions  are  computed  in  this 
fcishion.  The  result  is  an  algorithm  which- computes  the  field  due  to  n  particles 
at  n  locations  in  order  0(n)  operations. 

As  a  second  example,  the  evaluation  of  Legendre  expansions  can  be  accom¬ 
plished  by  application  of  a  particular  dense  upper-triangular  matrix  to  an  ar¬ 
bitrary  vector.  Away  from  the  main  diagonal,  the  elements  of  the  matrix  vary 
smoothly  cis  a  function  of  their  indices.  Alpert  and  Rokhlin  [2]  constructed  an 
algorithm  by  dividing  the  matrix  into  submatrices,  each  of  which  is  separated 
from  the  diagonal,  as  shown  in  Fig.  1.2.  This  division  ensures  that  each  subma¬ 
trix  has  smoothly  varying  elements  which  can  be  approximated  by  a  low-order 
polynomial.  The  error  of  the  approximation  decays  exponentially  in  the  degree 
of  the  polynomial,  so  in  practice  any  desired  degree  of  accuracy  may  be  obtained. 
The  algorithm  for  application  of  the  nxn  matrix  to  a  vector  then  requires  0(n) 
operations,  as  do  the  algorithms  of  [8],  [9],  [15]. 

Two  deficiencies  of  these  algorithms  are  apparent.  First,  though  they  share  a 


Figure  1.1:  The  field  due  to  particles  inside  circle  A  can  be  represented  by  a 
Laurent  expansion  about  the  center  of  A.  The  expansion  converges  rapidly  outside 
B,  which  has  radius  twice  that  of  A,  so  only  several  terms  of  the  expansion  are 
needed  to  evaluate  (to  high  precision)  the  effect  of  particles  inside  A  on  particles 
inside  C. 
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Figure  1.2:  Dense  matrix  connecting  coefficients  of  Legendre  and  Chebyshev  ex¬ 
pansions  is  divided  so  that  each  submatrix  can  be  approximated  to  high  accuracy 
by  a  low-order  polynomial.  This  combination  of  divide- and- conquer  and  efficient 
approximation  results  in  an  algorithm  to  apply  the  n  x  n  matrix  to  an  arbitrary 
vector  in  order  0(n)  operations. 

common  thread,  each  algorithm  is  custom-designed  for  its  own  problem  and  is  not 
generally  applicable  to  others.  Second,  in  the  case  when  the  matrix  is  invertible,  it 
is  not  clear  from  these  algorithms  how  to  efficiently  invert  the  matrix  or  apply  its 
inverse  to  a  vector.  These  deficiencies  have  recently  been  addressed  by  Beylkin, 
Coifman,  and  Rokhlin  in  [3], 


1.2  Wavelet  Bases 

1.2.1  Function  Space  Bases 

In  this  thesis  we  construct  a  clciss  of  orthonormal  bases  of  £^[0,  l]  (and  other 
spaces),  with  the  property  that  “smooth”  linear  operators  in  these  baises  corre¬ 
spond  to  matrices  which  are  sparse,  to  high  precision.  This  clciss  of  operators 
includes  those  mentioned  above,  as  well  as  operators  resulting  from  a  wide  va¬ 
riety  of  second-kind  integral  equations.  The  matrices  are  sparse;  furthermore, 
their  inverses  are  sparse,  to  high  precision.  As  a  result,  the  Schulz  method  [17], 
which  is  an  iterative,  quadratically  convergent  technique  for  matrix  inversion, 
becomes  highly  effective  and  produces  an  order  0(n  log^  n)  algorithm  for  solving 
a  second-kind  integral  equation. 
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Figure  1.3;  Multi-wavelet  bases  of  £^[0,1]  consist  of  functions  non-zero  on  in¬ 
tervals  of  various  scales.  The  basis  of  order  k  contains  k  functions  supported  on 
each  interval  of  one  scale,  which  does  not  overlap  other  intervals  of  that  scale. 
The  total  number  of  functions  on  a  single  scale  is  shown. 

A  matrix  is  sparse  to  high  precision  if  it  is  a  sparse  matrix  perturbed  by  a 
small  matrix,  i.e.,  if  most  of  the  elements  are  negligible.  In  actual  computations 
the  smallest  elements  are  discarded,  according  to  a  threshold  determined  by  a 
user-specified  accuracy. 

What  are  the  properties  of  a  basis  of  wavelets  that  lead  to  sparse  representa¬ 
tions?  The  baises  we  construct  are  orthonormal  and  consist  of  functions  that  are 
(generally) 

1.  orthogonal  to  low-order  polynomials,  and 

2.  non-zero  on  finite  intervals  of  various  lengths. 

A  basis  function  b  :  [0,  l]  — ♦  7^  is  orthogonal  to  low-order  polynomials  if  its  first 
k  moments  vanish. 


for  some  positive  integer  k.  In  our  construction,  ail  but  k  basis  functions  have  this 
property.  The  property  of  support  on  intervals  of  various  lengths  is  illustrated  in 
Fig.  1.3.  These  two  properties  ensure  that  a  function  that  is  smooth,  except  for 
a  finite  number  of  singularities,  will  have  a  negligible  projection  on  most  basis 
functions. 
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1.2.2  Vector  Space  Bases 

Numerical  solution  of  an  integral  equation,  as  mentioned  above,  can  be  accom¬ 
plished  by  projecting  onto  an  n-dimensional  subspace  of  and  solving  the  re¬ 
sulting  system  of  n  equations  in  n  unknowns.  Consider  the  second-kind  integral 
equation 

(I -)€)/  =  9, 

where  AC  is  an  integral  operator  given  by  the  formula 

(AC/)(i)=  I' KixN)f{t)dt, 

Jo 

and  suppose  that  {6i, . . . ,  6„}  is  a  orthonormal  basis  for  the  subspace.  Then  it  is 
easy  to  see  that  the  projection  is  obtained  by  computing  the  integrals 

f  f  K(x,t)  bi{x)  bj{t)  dx  dt,  i,;  =  l,...,n.  (1.1) 

Jo  Jo 

VVe  demonstrate  this  method  of  solution,  using  the  new  bases,  in  the  next  chapter. 
For  many  kernels,  however,  difficulties  may  arise  in  obtaining  the  integrals  (1.1). 

An  alternative  general  method  for  solution  of  integral  equations,  developed 
by  Nystrom,  approximates  the  integral  operator  by  a  quadrature,  resulting  in  a 

matrix  A  =  {wj  K{xi,  Xj)}i,j=\ . n  for  discretization  points  {xi, . . . ,  x„}  C  [0,  l| 

and  quadrature  weights  {t^i, . . .  The  matrix  A  is  dense,  co  by  direct  meth¬ 
ods  its  application  to  a  vector  requires  order  O(n^)  operations  and  its  inversion 
requires  O(n^)  operations. 

These  time  complexities  can  be  cut  considerably.  Complementary  to  the 
function  space  ba^es  sketched  above,  we  also  develop  a  vector  space  analogue,  for 
the  n-dimensional  space  of  functions  defined  on  {xi,...,Xn}.  The  vector  space 
bcises  have  properties  entirely  analogous  to  the  function  space  bases:  namely, 
all  but  k  bcisis  vectors  are  orthogonal  to  the  moment  vectors  {xi-', . . . ,  x„-’}  for 
j  =  0, . . . ,  ^  —  1,  and  the  basis  vectors  are  non-zero  on  various  “local”  subsets 
of  {xi,...,x„}.  In  solving  the  integral  equation  by  the  Nystrom  method,  the 
matrix  A  is  similarity-transformed  by  the  vector  basis  to  a  sparse  matrix  in 
order  O(nlogn)  operations.  In  addition,  its  inverse  is  sparse,  and  is  obtained  in 
order  0{n  log^  n)  operations  via  Schulz  method  (as  with  the  case).  Thus  the 
function  space  and  vector  space  bases  are  close  analogues,  and  both  lead  to  fast 
algorithms  for  solving  a  wide  cl«iss  of  integral  equations. 


1.3  Thesis  Organization 

In  Chapter  2  we  construct  the  function  space  bases  and  demonstrate  their  com¬ 
pleteness  and  convergence  properties.  We  prove  that  a  variety  of  integral  opera¬ 
tors  and  their  inverses,  when  expanded  in  these  bases,  are  represented  as  sparse 
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matrices.  We  present  analytical  expressions  for  the  bases  and  several  numerical 
examples  of  their  application  to  the  solution  of  second-kind  integral  equations. 

In  Chapter  3  we  construct  quaaratures  for  operator  kernels  with  certain  types 
of  diagonal  singularities.  The  quadratures  have  the  property  that  they  make 
only  endpoint  adjustments  to  the  trapezoidal  rule,  and  therefore  preserve  the 
smoothness  properties  of  the  kernel,  yet  they  achieve  high-order  convergence. 

In  Chapter  4  we  construct  the  vector  space  bases,  an  algorithm  for  their 
computation,  and  an  algorithm  to  transform  a  matrix  representing  an  integral 
operator  with  a  diagonally-singular  kernel  to  sparse  form  in  order  O(nlogn) 
operations.  We  employ  the  quadratures  of  Chapter  3  to  obtain  the  solution  of 
several  sample  integral  equations.  We  also  demonstrate  the  effectiveness  of  the 
new  bases  with  timing  results. 

Finally,  in  Chapter  5  we  mention  some  other  fast  algorithms  for  the  numerical 
solution  of  integral  equations,  we  discuss  extensions  to  the  methods  described  in 
the  preceding  chapters,  and  we  present  further  applications  for  the  new  bases. 


Chapter  2 

Function  Space  Bases 


Families  of  functions  ha,bi 

a,b  £  TZ,  a  ^  0, 

derived  from  a  single  function  h  by  dilation  and  translation,  which  form  a  beisis 
for  are  known  as  wavelets  (Grossman  and  Morlet  [10]).  In  recent  years, 

these  families  have  received  study  by  many  authors,  resulting  in  constructions 
with  a  variety  of  properties.  Meyer  [13]  constructed  orthonormal  wavelets  for 
which  h  €  C°^{1Z).  Daubechies  [5]  constructed  compactly  supported  wavelets 
with  h  €  C'‘{TZ)  for  arbitrary  k,  and  [5]  gives  an  overview  and  synthesis  of  the 
field. 

In  this  chapter  we  construct  a  somewhat  different  type  of  basis  for  C^{TZ)  that 
can  be  readily  revised  to  a  basis  for  £^[0, 1].  Each  basis,  which  we  call  a  multi¬ 
wavelet  basis,  is  comprised  of  dilates  and  translates  of  a  finite  set  of  functions 
/ii, . . . ,  /iji;.  In  particular,  our  bcises  consist  of  orthonormal  functions 


haM  =  |a| 


hl^{x)  =  2'"/^  h,i2^x  -  n),  j  =  h...,k;m,nEZ,  (2.1) 

where  hi,...,hk  are  piecewise  polynomial,  are  supported  on  the  interval  [0,1], 
and  have  vanishing  moments, 

f  hj{x)  x' dx  =  0,  1  =  0,1,...,^  —  !.  (2.2) 

Jo 

The  properties  of  compact  support  and  vanishing  moments  lead  to  bases  in  which 
a  variety  of  integral  operators  are  represented  by  sparse  matrices.  In  particular, 
an  integral  operator  whose  kernel  is  non-oscillatory  and  analytic  e.xcept  along  a 
finite  set  of  curves,  when  expanded  in  one  of  these  bases,  is  sparse. 

In  the  following  sections,  we  construct  multi-wavelet  bases  in  one  and  several 
dimensions,  prove  their  completeness,  convergence,  and  sparsity  properties,  and 
present  several  examples  to  illustrate  their  numerical  usefulness. 
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2.1  Multi- Wavelet  Bases 

2.1.1  The  One-Dimensional  Construction 

We  first  restrict  our  attention  to  the  finite  interval  [0, 1]  C  ^  and  we  construct  a 
basis  for  £^[0,1].  We  employ  the  multi-resolution  analysis  framework  developed 
by  Mallat  [12]  and  Meyer  [14],  and  discussed  at  length  by  Daubechies  [5].  We 
suppose  that  A:  is  a  positive  integer  and  for  m  =  0, 1, 2, . . .  we  define  a  space  S^ 
of  piecewise  polynomial  functions, 

5^  =  {/  :  the  restriction  of  /  to  the  interval  (2~”*n, 2“"*(n  -I-  1))  is  (2.3) 
a  polynomial  of  degree  less  than  A:,  for  n  =  0, . . . ,  2"*  —  1, 
and  /  vanishes  elsewhere}. 

It  is  apparent  that  the  space  has  dimension  2”*  A;  and 

So"c5{=C---C5^C---. 

For  m  =  0, 1, 2, ...  we  define  the  2"*  A:-dimensional  space  R'^  to  be  the  orthogonal 
complement  of  5^  in  S^^i, 

=  si*,.  .Rixsi, 

so  we  have  the  decomposition 

©  /I*  0  •  • .  ©  (2.4) 

Suppose  that  functions  h\, . . .  ,hk  :  TZ  TZ.  form  an  orthogonal  basis  for  Rq. 
By  the  orthogonality  of  Rq  to  Sq,  the  first  k  moments  of  hi,. . .  ,hk  vanish, 

f  hj{x)  x' dx  =  0,  i  =  0, 1, . . . ,  A:  —  1. 

Jo 

The  2l'-dimensional  space  iZj  is  spanned  by  the  2A:  orthogonal  functions  /ii(2x), 
. . . ,  hk{2x),  h\{2x  —  1), . . . ,  hk{2x  —  1),  of  which  k  are  supported  on  the  inter¬ 
val  [0,^]  and  k  on  [^,1].  In  general,  the  space  is  spanned  by  2”* A;  functions 
obtained  from  hi,...,hic  by  translation  and  dilation.  There  is  some  freedom  in 
choosing  the  functions  hi,...,hk  within  the  constraint  that  they  be  orthogo¬ 
nal;  by  requiring  normality  and  additional  vanishing  moments,  we  specify  them 
uniquely,  up  to  sign. 

In  preparation  for  the  definition  of  hi,...,hk,  we  construct  the  k  functions 
fii  ■  ■  ■  1  fk  '•  TZ  —*  TZ,  supported  on  the  interval  [—1, 1],  with  the  following  prop¬ 
erties: 

1.  The  restriction  of  fi  to  the  interval  (0, 1)  is  a  polynomial  of  degree  A-  —  1. 
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2.  The  function  fi  is  extended  to  the  interval  (—1, 0)  as  an  even  or  odd  function 
according  to  the  parity  of  i  +  A:  —  1. 

3.  The  functions  satisfy  the  following  orthogonality  and  normality 

conditions: 

J  ^  /«(^)  fji^)  =  {fit  fj)  ~  ^ij  ~  !»•••) 

4.  The  function  fj  has  vanishing  moments, 

J  fj{x)x'dx  =  0,  i  =  0,1,. . .  ,j  +  k  -  2. 

Properties  1  and  2  imply  that  there  are  polynomial  coefficients  that  deter¬ 
mine  the  functions  fi,...,fki  while  properties  3  and  4  provide  k^  (non-trivial) 
constraints.  It  turns  out  that  the  equations  uncouple  to  give  k  nonsingular  lin¬ 
ear  systems  that  may  be  solved  to  obtain  the  coefficients,  yielding  the  functions 
uniquely  (up  to  sign).  Rather  than  prove  that  these  systems  are  nonsingular, 
however,  we  now  determine  fi,...,fk  constructively. 

We  start  with  2k  functions  which  span  the  space  of  functions  that  are  poly¬ 
nomials  of  degree  less  than  k  on  the  interval  (0,1)  and  on  (—1,0),  then  or- 
thogonalize  k  of  them,  first  to  the  functions  1,  x, . . . ,  x*'"^  then  to  the  functions 
X*,  X*'*’* , . . . ,  x^*~^  and  finally  among  themselves.  We  define  fl  t>y  the 

formula 

[  xJ"S  X  €  (0,1), 

/](a:)  =  S  xe(-l,0), 

[  0,  otherwise, 

and  note  that  the  2k  functions  1 ,  x, . . . ,  x*”‘ ,  are  linearly  indepen¬ 

dent,  hence  span  the  space  of  functions  that  are  polynomials  of  degree  less  than 
k  on  (0, 1)  and  on  (  —  1,0). 

1.  By  the  Gram-Schmidt  process  we  orthogonalize  fj  with  respect  to  1,  x, . . . , 
x*“',  to  obtain  fJ,  for  j  =  1,...,/;.  This  orthogonality  is  preserved  by  the 
remaining  orthogonalizations,  which  only  produce  linear  combinations  of 
the  fJ. 

2.  The  next  sequence  of  steps  yields  k  —  I  functions  orthogonal  to  x*',  of 

which  k  —  2  functions  are  orthogonal  to  x*'*’*,  and  so  forth,  down  to  1 
function  which  is  orthogonal  to  First,  if  at  least  one  of  fJ  is  not 

orthogonal  to  x^,  we  reorder  the  functions  so  that  it  appears  first,  {f^,x'‘)  ^ 
0.  We  then  define  ff  =  fJ  —  aj  ■  where  Oj  is  chosen  so  (/?,x^)  =  0 
for  j  =  2, . . . ,  k,  achieving  the  desired  orthogonality  to  x*'.  Similarly,  we 
orthogonalize  to  each  in  turn,  to  obtain  /j,  •  •  •  >  fk^^ 

such  that  (/j^^x*)  =  0  for  i  <  j  -|-  Jt  —  2. 
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3.  Finally,  we  do  Gram-Schmidt  orthogonalization  on  •  •  •  >/?>  in 

that  order,  and  normalize  to  obtain  /*, /k-i, .  fi¬ 
ll  is  readily  seen  that  the  fj  satisfy  properties  1-4  of  the  previous  paragraph. 
Defining  hi,.  ..,hk  :  TZ  —*  Hhy  the  formula 

we  obtain  the  equality 

=  linear  span  {h,  :  i  =  .  ,k}, 

and,  generally, 

Rm  =  linear  span  :  hl„,{x)  =  2’"/^  hj{2’^x  -  n),  . 

j  =  l,...,k‘,  n  =  0,...,2--l}.  ^  •  ’ 

We  will  show  next  that  dilates  and  translates  of  the  piecewise  polynomial  func¬ 
tions  hi,...,hk  form  an  orthonormal  basis  for  C^{'JZ).  Furthermore,  a  subset 
of  these  dilates  and  translates,  combined  with  a  basis  for  Sq,  forms  a  basis  for 

£»:o,  1]. 

2.1.2  Completeness  of  One-Dimensional  Construction 

We  define  the  space  S*'  to  be  the  union  of  the  given  by  the  formula 

s‘  =  U  SL  (2.6) 

m=0 


and  observe  that  S*'  =  £^[0, 1].  In  particular,  contains  the  Haar  basis  for 
£^[0, 1],  consisting  of  functions  piecewise  constant  on  each  of  the  subintervals 
(2“”*n,  2"”*(n  -I-  1)).  Here  the  closure  S*'  is  defined  with  respect  to  the  £^-norm, 

ll/ll  =  (/./>■'^ 

where  the  inner  product  {f,g)  is  defined  by  the  formula 

{f,9)  =  [  f{x)g{x)dx. 

Jo 

We  let  ux, . . .  denote  an  orthonormal  baisis  for  Sq-,  in  view  of  Eqs.  (2.4),  (2.5), 
and  (2.6),  the  orthonormal  set 

{ u  j  I  j  —  1  ^  At  ^ 

u  {^i,m  :  j  =  I,. ..  ,k;  m  =  0,1,2,...;  n  =  0, . . . ,  2”*  -  1} 
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spans  £^[0,  l];  we  refer  to  it  eis  the  multi-wavelet  basis  of  order  k  for  £*[0,  l]. 

Now  we  construct  a  basis  for  £^(72.)  by  defining,  for  m  €.  Z,  the  space  by 
the  formula 

=  {/  :  the  restriction  of  /  to  the  interval  (2“"‘n,  2“”*(n  +  1))  is 
a  polynomial  of  degree  less  than  A:,  for  n  G  2} 

and  observing  that  the  space  is  spanned  by  the  orthonormal  set 

=  2"'^  A,(2”x  -  n),  j  =  1, . . . ,  n  6  2}. 


Thus  which  is  contained  in  orthonormal  basis 

{^j.m  :  i  =  1,...,A:;  m,n  6  2}. 


2.1.3  Construction  in  Multiple  Dimensions 

The  construction  of  our  bases  for  £^[0, 1]  and  £^(72.)  can  be  extended  to  certain 
other  function  spaces,  including  £^[0,6]*^  and  £^(7^“^),  for  any  positive  integer  d. 
We  now  outline  this  extension  by  giving  the  basis  for  £^[0, 1]*,  which  is  illustrative 
of  the  construction  for  any  finite-dimensional  space.  We  define  the  space  S'^  by 
the  formula 

S‘’“  =  S‘xS‘,  m  =  0,l,2 . 

where  S’^  is  defined  by  Eq.  (2.3).  We  further  define  R}^  to  be  the  orthogonal 
complement  of  S^  in  5^’^i , 


fc.2 


Jk  2  • 

Then  R^'  is  the  space  spanned  by  the  orthonormal  basis 


{ui{x)hj{y),  hi{x)uj{y),  hi{x)hj{y)  :  ij  =  l,...,k}. 

Among  these  basis  elements  each  element  u{x,y)  has  no  projection  on  low- 
order  polynomials, 

nv{x,y)x' t/-'  dxdy  =  0,  i,;'  =  0, 1, . . . ,  A:  -  1. 

.1 

The  space  is  spanned  by  dilations  and  translations  of  the  v{x,  y)  and  the  basis 
of  £^[0, 1]^  consists  of  these  functions  and  the  low-order  polynomials  {ui{x)uj{y)  : 
i,j  =  1,...,*:}. 
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2.2  Convergence  of  the  Multi- Wavelet  Bases 


For  a  function  /  €  ^^[0,  l],  a  positive  integer  k,  and  m  =  0, 1, 2 . . we  define  the 
orthogonal  projection  Q^f  of  /  onto  by  the  formula 

j,n 

where  is  an  orthonormal  basis  for  S^.  The  projections  Q'^f  converge  (in 

the  mean)  to  /  as  m  — >  oo.  If  the  function  /  is  several  times  differentiable,  we 
can  bound  the  error,  as  established  by  the  following  lemma. 

Lemma  2.1  Suppose  that  the  function  f  :  [0, 1]  — >  72.  is  k  times  continuously 
differentiable,  f  €  C*'[0, 1].  Then  Q'^f  approximates  f  with  mean  error  bounded 
as  follows: 

Kif  -  /II  <  2-"‘  sup  |/l‘>(x)|  (2.7) 

4**:!  igfo.i) 

Proof  We  divide  the  interval  [0, 1]  into  subintervals  on  which  Q^f  is  a  poly¬ 
nomial;  the  restriction  of  Q^f  to  one  such  subinterval  Im,n  is  the  polynomial  of 
degree  less  than  k  that  approximates  /  with  minimum  mean  error.  We  then  use 
the  maximum  error  estimate  for  the  polynomial  which  interpolates  /  at  Cheby- 
shev  nodes  of  order  k  on  Im,n- 

We  define  /m.n  =  (2“"‘n,2'"”(n 1)]  for  n  =  0, 1, . . .  ,2*"  -  1,  and  obtain 


WQif-ff 


=  I' [(Qtf)M  - /mY 

=  1]/,  [(<3m/)(^)  - 

n  *m,n 

<  Ef  [(Cl„/)(x)-/(x))"<ix 

n  ^  lm,n 


< 


sup  l/<*>(x)| 


2 

dx 


Ol-mA: 

-rrrr  sup  |/‘*'^(x)l 
4*«r.'  r€[0,lj 


and  by  taking  square  roots  we  have  bound  (2.7).  Here  C^„f  denotes  the  poly¬ 
nomial  of  degree  k  which  agrees  with  /  at  the  Chebyshev  nodes  of  order  k  on 
fm.n,  and  we  have  used  the  well-known  maximum  error  bound  for  Chebyshev 
interpolation  (see,  e.g.,  [4]).  □ 

The  error  of  the  approximation  Q'^f  of  /  therefore  decays  like  2"'"*'  and, 
since  S!^  has  a  basis  of  2"'k  elements,  we  have  convergence  of  order  k.  For  the 
generalization  to  d  dimensions,  it  is  easily  seen  that  the  rate  of  convergence  is  of 
order  k/d. 
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2.3  Second-Kind  Integral  Equations 

A  linear  Fredholm  integral  equation  of  the  second  kind  may  be  written  in  the 
form 

f{x)-  [  K{x,t)f{t)dt  =  g{x),  (2.8) 

where  the  function  K  is  the  kernel,^  is  the  right-hand-side,  and  /  is  the  unknown. 
For  notational  simplicity,  in  this  chapter  we  restrict  our  attention  to  the  interval 
[a,  6]  =  [0, 1].  We  use  the  symbol  JC  to  denote  the  integral  operator  of  Eq.  (2.8), 
given  by  the  formula 

{Kf){x)=  [' K{x,t)f{t)dt, 

Jo 

for  all  /  6  £^[0,1]  and  x  €  [0,1].  Suppose  that  {bi,b2,...}  is  an  orthonormal 
basis  for  £^[0,  Ij;  the  expansion  of  K  in  this  basis  is  given  by  the  formula 

Kix,t)  =  '£f2K,,bi{x)bj{t),  (2.9) 

t=i  j=i 

where  the  coefficient  Kij  is  given  by  the  expression 

Kij  =  f  f  K{x,t)  bi{x)  bj{t)  dx  dt,  i,j=:l,2, —  (2.10) 

Jo  Jo 

Similarly,  the  functions  /  and  g  have  expansions 

/(®)  =  £  /.  ^(3;).  9{x)  =  f^9i  biix), 

i=l  i=l 

where  the  coefficients  fi  and  gi  are  given  by  the  formulae 

/.  =  /  f{x)bi{x)dx,  9i=  [  g{x)bi{x)dx,  i  =  l,2,.... 

Jo  Jo 

The  integral  equation  (2.8)  then  becomes  the  system  of  equations 

fi~fll<i:fo=9i,  i  =  l,2,.... 

j=i 

The  expansion  for  K  may  be  truncated  at  a  finite  number  of  terms,  yielding  the 
integral  operator  R  defined  by  the  formula 

(fl/)(i)  =  f' EE  (K-i  I’M  ijW)  m 

t=0  j=0 


/€£'[0,1],  x€[0,l]. 
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which  approximates  tC.  Integral  equation  (2.8)  is  thereby  approximated  by  the 
system 

i  =  (2.11) 

a  system  of  n  equations  in  n  unknowns.  Eqs.  (2.11)  may  be  solved  numerically 
to  yield  an  approximate  solution  to  Eq.  (2.8),  given  by  the  expression 

/n(i)  =  E  /. 

t=l 

How  large  is  the  error  en  =  /  —  /r  of  the  approximate  solution?  We  follow 
the  derivation  by  Delves  and  Mohamed  in  [6].  Rewriting  Eqs.  (2.8)  and  (2.11)  in 
terms  of  operators  AC  and  we  have 

=  9 

{I-R)fR  =  g, 

and  combining  the  latter  equations  yields 

(/  -  K)eR  =  (AC  -  R)fR. 

Provided  that  (/  —  AC)“*  exists,  we  obtain  the  error  bound 

(2-12) 

The  error  depends,  therefore,  on  the  conditioning  of  the  original  integral  equation, 
cis  is  apparent  from  the  term  ||(/  —  AC)~‘||,  and  on  the  fidelity  of  the  finite¬ 
dimensional  operator  R  to  the  integral  operator  AC. 


2.4  Sparse  Representation  of  Integral  Opera¬ 
tors  and  Their  Inverses 

2.4.1  Representation  in  Multi- Wavelet  Bases 

We  consider  integral  operators  AC  with  kernels  that  are  analytic,  except  at  x  = 
t,  where  they  are  singular.  In  particular,  we  analyze  singularities  of  the  form 
log  |x  —  <1  or  the  form  |x  —  t|“,  with  0  <  |q|  <  1.  An  operator  with  such  a  kernel 
K,  expanded  in  one  of  the  multi-wavelet  bases  defined  above,  is  represented  as 
a  sparse  matrix.  This  sparseness  is  due  to  the  smoothness  of  K  on  rectangles 
separated  from  the  “diagonal”. 
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Figure  2.1:  Rectangular  regions  (just)  separated  from  the  diagonal. 

Definition  2.2  We  say  that  a  rectangular  region  oriented  parallel  to  the  coor¬ 
dinate  axes  X,  t  is  separated  from  the  diagonal  if  its  distance  in  the  horizontal  or 
vertical  direction  from  the  line  x  =  <  is  at  least  the  length  of  its  longer  side.  In 
symbols,  a  region  [x,x  -I-  a]  x  [t,  t  -f  6]  C  TV  is  separated  fi’om  the  diagonal  if 
a  -}-  max{a,  b}  <  t  —  x  or  6  -f-  max{a,  b}  <  x  —  t. 

This  definition  is  illustrated  in  Fig.  2.1. 

Suppose  that  k  is  a.  positive  integer  and  that  {6i,  62,  •  •  is  the  multi-wavelet 
basis  for  £^[0, 1]  of  order  k,  defined  in  Section  2.1.  We  let  Ij  denote  the  interval 
of  support  of  6j,  and  we  <issume  that  the  sequence  of  beisis  functions  61,  621  •  •  •  is 
ordered  so  that  /i,/2, ...  have  non-increasing  lengths.  For  large  n,  the  matrix 
. n  is  sparse,  to  high  precision,  as  is  proved  in  the  following  proposi¬ 
tions. 

Lemma  2.3  Suppose  that  the  function  K  :  [0, 1]  x  [0, 1]  — >  7^  is  given  by  the 
formula  K{x,t)  =  log  |x  —  <|.  The  expansion  of  K  in  the  multi-wavelet  basis  of 
order  k  (Eq.  2.9)  has  coefficients  Kij  which  satisfy  the  bound 

|A'.,I  <  (2.13) 

whenever  the  rectangular  region  L  x  Ij  is  separated  from  the  diagonal. 

Proof.  Suppose  that  the  intervals  U  and  Ij  are  given  by  the  expressions 
/,  =  -f  a]  and  Ij  =  [to,  to  +  6];  without  loss  of  generality  we  assume  (one 
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of  two  equivalent  cases)  that  b  +  max{a,6}  <  xo  ~  ^o-  It  is  immediate  from  this 
inequality  that 

xq  +  af2  —  X 
Xq  +  a/2  —  t 

for  (x,/)  €  It  X  Ij- 

We  use  the  Taylor  expansion  for  the  natural  logarithm  about  c  >  0, 

Iog(c  -i-y)  =  log(c)  +  (y/c)  ~  {yfcfi2  +  (y/cf/3  -  (ylc^fi  +  •  •  • , 


for  |j/|  <  c.  We  now  let  c  =  xq  +  a/2  —  t  and  y  =  x  —  xq  ~  a/2  and  for 
(x,<)  &  li  X  Ij  we  obtain  the  formula 


log  |x  -  tl  =  log(xo  +  a/2  -  t)  -  ^ 

V^o  +  a/2- 


(2.15) 


We  now  apply  Eqs.  (2.10),  (2.15),  (2.2),  and  (2.14),  each  in  turn,  to  obtain 


=  /  /  K{x,t)b,{x)bj{i)dx  dt 

|./to  Jxo 

fta+b  I  /To+a 

<  y  log  |x  —  t|  6i(x)  dx  \bj{t)\dt 

Jfto+b  I  fxo+a  n 

[  \L  K»+2-‘> 

fto+b  fxo+a  \  »  /IN'” 

^  i  L  i£u) 

fto+b  fXo+a  1 

-  C  'JiC'  UT  ’ 


1  rto+b  /  /  rxo+a 

-  2k  ■  3*-l  L  yUro 

<  — <  _1 _ 

-  2it-3*-*  -  Sit -a*-!’ 


as  was  to  be  proved.  □ 


Lemma  2.4  Suppose  that  the  function  L  :  D  x  D  —*  C  is  given  by  the  formula 
L{z,w)  —  f{z,w)[og\2  —  w\,  where  D  is  the  closed  disk  of  radius  |  centered  at 
2  =  1  and  f  is  analytic  in  a  domain  containing  D  x  D  C.  C^.  Suppose  further 
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that  the  function  K  is  the  restriction  of  L  to  [0,1]  x  [0,1].  The  expansion  of  K 
in  the  multi-wavelet  basis  of  order  k  has  coefficients  Kij  which  satisfy  the  bound 


whenever  the  rectangular  region  li  x  Ij  is  separated  from  the  diagonal. 


(2.16) 


Proof.  We  combine  the  method  of  proof  used  in  Lemma  2.3  with  the  formula 
for  the  derivative  of  a  product. 


dx”^  ^  V  / 

By  the  Cauchy  integral  formula  we  obtain 


t)  ’■  log  |x  —  t| 


dt' 


<  r!  sup  1/(2,  ti;) 

z.xo^dD 


for  (x,<)  €  [0,  l]  X  [0, 1].  For  the  logarithm,  differentiation  yields  the  formula 
a”‘-Mog|x-t|  (-l)"*-’-‘(m-r-l)! 

~  ^3.  _  lyn-T  ’ 

for  r  <  m.  Combining  these  expressions,  we  obtain 


d”'I<{x,t) 

VI 

(rn^ 

f{x,t) 

log  |x  —  t| 

dx”' 

Kn 

dx^ 

5x"‘-’’ 

<  sup  |/(2,u;)|  +  m!  [log  |x  -  t|| 

^  C  2  +  logm\ 


<  S/  •  (  m! 


|x  -  t\^ 


(2.17) 


for  |x  —  <  1  and  m  >  1,  where  5/  =  l/(^i«^)l- 

Suppose  that  the  intervals  L  and  Ij  are  given  by  the  expressions  /,  =  [xq,  xq  4- 
a]  and  Ij  =  [to,  to  +  &];  we  assume  without  loss  of  generality  that  6  +  max{a,  b)  < 
xo  —  to-  It  follows  directly  from  this  inequality  that 


JQ  +  q/2  —  X  ^  1^ 


(2.18) 


1  Xo  +  a/2  -  1 1  -  3  '  ^ 

for  (x,t)  e  /,  X  Ij.  We  now  apply  Eqs.  (2.10),  (2.2),  (2.17),  and  (2.18),  to  obtain 
I  rto+b  rxo+a  I 


rto+0  AXo+a 

=  /  /  K{x,t)bi{x)bj{t)dx 

Jto  Jtq 


18 


CHAPTER  2.  FUNCTION  SPACE  BASES 


< 

< 

< 

< 

< 

< 


rto+k  I  rio+a  »  (jQ  +  a/2  -  x)"*  d”'K{xo  +  q/2,0  ^ 

Jto  \Jxo  “o 

rto+b  fx 

Jto  Jxa 


\bj{t)\dt 


r<o+fc  yro+a  ^  Xq  +  af2  —  X 

Lo  ^0  +  al2-t 

r<o+^  r*o+* 


5/  (2  +  logm)  \bi{x)\ dx  \bj{t)\ dt 


rto+  /•XO+  g  /1\  ^ 

•/«0  -/xn  TO=fc 


which  was  to  be  proved.  0 

The  proofs  of  the  following  two  lemmiis  closely  resemble  those  of  Lemma  2.3 
and  Lemma  2.4,  and  are  omitted. 

Lemma  2.5  Suppose  that  the  function  K  :  [0,  l]  x  [0,  \\  —*  'R-  is  given  by  the 
formula  I\{x,t)  =  [x  -  <|®  with  0  <  |al  <  1.  Then  the  expansion  coefficient  Ki, 
of  the  function  K  in  the  multi-wavelet  basis  of  order  k  satisfies  the  bound 

(219) 

whenever  the  rectangular  region  /,  x  Ij  is  separated  from  the  diagonal. 


Lemma  2.6  Suppose  that  the  function  L  :  D  x  D  —*  C  is  given  by  the  formula 
L{z,  w)  =  f{z,  u;)|2-u;|°',  with  0  <  |a|  <  I,  where  D  is  the  closed  disk  of  radius  | 
centered  ai  r  =  ^  and  f  is  analytic  in  a  domain  containing  D  x  D  C  C^-  Suppose 
further  that  the  function  K  is  the  restriction  of  L  to  [0, 1]  x  [0,  l].  The  expansion 
of  K  in  the  multi-wavelet  basis  of  order  k  has  coefficients  hij  which  satisfy  the 
bound 


\i<,A  < 


1 

3*-‘ 


sup  \f{z,w)\, 

z.w^dD 


(2.20) 


whenever  the  rectangular  region  7,  x  Ij  is  separated  from  the  diagonal. 


The  four  preceding  lemmas  show  that  for  a  smooth  kernel  /v  with  logarithm 
or  power  singularity  at  x  =  t,  the  order  k  of  the  multi-wavelet  basis  in  which 
K  is  expanded  may  be  chosen  large  enough  that  the  expansion  coefficient  h,j  is 
negligible,  provided  /,  x  Ij  is  separated  from  the  diagonal.  X  similar  statement 
can  be  proven  for  any  kernel  of  the  form  /\'(x,  t)  =  (^(x,  f)s(|x  — fD-t-t/iix,  <),  where 
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are  entire  analytic  functions  of  two  variables  and  s  is  an  analytic  function 
except  at  the  origin  (where  it  has  a  singularity),  provided  that  s  is  integrable. 
We  do  not  prove  this  statement  here. 

The  next  lemma  establishes  the  fact  that,  asymptotically,  most  regions  7,  x  Ij 
are  separated  from  the  diagonal. 

Lemma  2.7  Suppose  that  /i, . . . ,  7„  are  the  (non-increasing)  intervals  of  support 
of  the  first  n  functions  of  the  multi-wavelet  basis  of  order  k.  Of  the  rectangular 
regions  L  x  Ij,  we  denote  the  number  separated  from  the  diagonal  by  S{n)  and 
the  number  “near”  the  diagonal  by  N{n)  =  —  S(n).  Then  N{n)  grows  as 

O(nlogn);  in  particular,  for  n  =  2‘k  with  I  >  0,  we  have  the  formula 

N{n)  =  6nlk  -  ISnJfc  -  6lk^  +  ISk^.  (2.21) 

Proof.  The  restriction  that  n  =  2^k  ensures  that  the  first  n  basis  functions 
consist  of  those  functions  whose  intervals  of  support  have  length  at  least  2^~L 
We  define  5i(p)  to  be  the  number  of  pairs  {i,j)  such  that  the  rectangular  region 
/,  X  /j  is  separated  from  the  diagonal  and  |/ij  =  |/j|  =  2“^,  and  we  observe  that 

^i{p)  —  (2'*  —  1)(2^  —  2)  k^  for  p  =  0,1,2, _  We  further  define  S2(p,q)  to  be 

the  number  of  pairs  {i,j)  such  that  I,  x  Ij  is  separated  from  the  diagonal  and 
17, 1  =  2“’’,  |7j|  =  2*’,  and  we  observe  that  S2{p,q)  =  5i(min{p, 9})2l'’~’l  for 
p,  9  =  0, 1, 2, . . ..  Finally,  we  combine  these  formulae  to  obtain 

•5(«)  =  12  i^2{p,q)  +  S2{q,p)) 

p=0  \  <?=p+l 

=  £  S,(p)  (l  +  2(2'-'' -  2)) 

p=0 
(-1 

(2P-l)(2P-2)^'2(2'-'’+*  -3) 

p=0 

=  (4'  -  6  •  2‘l  +  15  •  2'  +  6/  -  16)  k^ 

=  —  6nlk -h  15nk -h  6lk^  —  16k^, 

from  which  Eq.  (2.21)  follows  directly.  The  assertion  that  the  general  growth  of 
N(n)  is  O(nlogn)  follows  from  Eq.  (2.21)  and  that  fact  that  is  a  monotonic 
function  of  n.  □ 

2.4.2  Products  of  Integral  Operators 

The  previous  subsection  established  the  fact  that  a  wide  class  of  integral  op¬ 
erators,  when  expanded  in  multi-wavelet  coordinates,  are  represented  to  high 
accuracy  m  sparse  matrices.  It  readily  follows  that  a  product  of  such  integral 
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operators  can  be  similarly  represented.  For  if  we  define  integral  operators  }Ci,K.2 
by  the  formulate 

(/Ci/)(x)  =  r  h\{x,t)f{t)dt 

Jo 

(^2/)(a:)  =  r  K2{x,t)f{t)dt, 

Jo 

then  the  product  operator  /C3  =  KL^Ki  is  given  by  the  formula 

{^2K.xf){x)  =  C  f  I<2{x,y)K,{y:l)f{t)dtdy 
Jo  Jo 

-  Jo  {Jo  /(O 

=  r  K^ix,t)f{t)dt, 

Jo 

where  the  kernel  K3  of  the  product  has  the  form 

Ksix^t)  =  f  K2{x,y)  Ki{y,t)dy. 

Jo 

If  kernels  Ki  and  K2  are  analytic  except  along  the  diagonal  x  =  t,  where  they 
have  integrable  singularities,  then  the  same  is  true  of  the  product  kernel  K3.  As  a 
result,  the  product  operator  K3  also  has  a  sparse  representation  in  a  multi- wavelet 
basis. 

2.4.3  Schulz  Method  of  Matrix  Inversion 

Schulz’s  method  [17]  is  an  iterative,  quadratically  convergent  algorithm  for  com¬ 
puting  the  inverse  of  a  linear  operator.  Its  performance  is  characterized  by  the 
following  lemma. 

Lemma  2.8  Suppose  that  A  is  an  invertible  linear  operator,  Xq  is  the  operator 
given  by  Xq  =  A^ l\\A^ A\\,  and  for  m  =  0, 1,2, . . .  the  operator  is  defined 

by  the  recursion 

Xffi^l  2Xfn  XfjiA-Xjn* 

Then  Xm+i  satisfies  the  formula 

1  -X^+iA  =  {I-X,nA)\  (2.22) 

Furthermore,  Xjn  —*  A~^  as  m  00  and  for  any  e  >  0  we  have 

||/-X,„/ll|<e  provided  m  >  2log2  «(A) -I- logj  log(l/e),  (2.23) 

where  k(A)  =  ||/l||  •  ||A~'||  is  the  condition  number  of  A  and  the  norm  ||.4|1  = 
(largest  eigenvalue  of  A^A)^^'^. 
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Proof.  Eq.  (2.22)  is  obtained  directly  from  the  definition  of  Xm+i-  Bound 
(2.23)  is  equally  straightforward.  Noting  that  A^ A  is  symmetric  positive-definite 
and  letting  Aq  denote  the  smallest  and  Ai  the  largest  eigenvalue  of  A^A  we  have 


\\I-XoA\\ 


A»A 

HAWAII 

1  —  Ao/ Ai 

1  -/c(A)"^ 


(2.24) 


From  Eq.  (2.22)  we  obtain  I  —  XmA  =  (/  —  XoA)^”,  which  in  combination  with 
Eq.  (2.24)  and  simple  manipulation  yields  bound  (2.23).  □ 

The  Schulz  method  is  a  notably  simple  scheme  for  matrix  inversion  and  its 
convergence  is  extremely  rapid.  It  is  rarely  used,  however,  because  it  involves 
matrix-matrix  multiplications  on  each  iteration;  for  most  problem  formulations, 
this  process  requires  order  O(n^)  operations  for  an  n  x  n  matrix.  As  we  have 
seen  above,  on  the  other  hand,  a  discretized  integral  operator  A  represented  in 
the  basis  of  Section  2.1  has  only  order  O(nlogn)  elements  (to  finite  precision). 
In  addition,  A^ A  and  (A^A)'"  are  similarly  sparse.  This  property  enables  us  to 
employ  the  Schulz  algorithm  to  compute  A~^  in  order  0(n  log^  n)  operations. 


2.5  Numerical  Examples 

2.5.1  Basis  Functions 

In  this  section  we  give  numerical  expressions  for  the  multi-wavelet  functions 
foi  fi,  •  •  •  1  fk-i  and  show  their  graphs  for  several  values  of  k.  Table  2.1  con¬ 
tains,  for  small  the  polynomials  which  represent  the  fi  on  the  interval  (0, 1), 
together  with  the  reflection  formula  to  extend  the  functions  to  (  —  1,1),  which  is 
their  interval  of  support.  Fig.  2.2  shows  the  graphs  of  the  functions  for  A:  =  4 
and  k  =  5. 

2.5.2  Integral  Operators  and  Their  Inverses 

We  compute  the  expansion  in  multi-wavelet  bases  of  the  integral  operator  K 
defin^’d  by  the  formula 


(AC/)(x)=  f\og\x-t\mdt, 
Jo 


(2.25) 


T  =  {  A', y 


which  yields  the  matrix 


(2.26) 
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Table  2.1:  Expressions  for  the  orthonormal,  vanishing-moment  functions 
fii‘--ifk,  for  various  k,  for  argument  x  in  the  interval  (0,1).  The  function 
fi  is  extended  to  the  interval  (  —  1,1)  as  an  odd  or  even  function,  according  to  the 
formula  fi{x)  =  (— 1)*'''*'“^/,(— x)  for  x  €  (—1,0),  and  is  zero  outside  (—1,1). 
The  functions,  given  here  for  fc  =  1, . . .  ,5,  are  tabulated  for  larger  values  of  k  in 
the  chapter  appendix. 


/l(x)  = 

A:  =  1 

J\ 

flix)  = 

k  =  2 

\l\  (-l+2x) 

Mx)  = 

(-2  +  3X) 

Mx)  = 

k  =  3 

(1  -24x  +  30i2) 

f2{x)  = 

\yj\  (3  -  16x -h  15x2) 

fsix)  = 

(4-15i  +  12i") 

fi{x)  = 

A:  =  4 

yg  (1 -f4x -30x2 -1-28x3) 

f2{x)  = 

^  (-4  H-  105x  -  300x2  210x3) 

/3(X)  = 

2\/il  (-5  +  48x  -  105x2  -h  64x3) 

(-16  -h  105x  -  192x2  -t-  105x3) 

f4{x)  = 

flix]  = 

A:  =  5 

yj  ( 1  -1-  30x  -1-  210x2  _  g40a;3  ^  ) 

f2ix)  = 

2\/5  (-5  -  144x  -i-  115.5x2  _  2240x3  -|-  1260x‘') 

hix)  = 

(22  -  735x  -I-  3504x2  _  5450x3  4.  2700x‘‘) 

U{x)  = 

8\/i  (35  -  512x  -h  1890x2  _  2560x3  4.  11553.4) 

fsix)  = 

2\/t£  (32  -315x -1-960x2  -  1155x3  -I- 480x'') 

where 

Kij  =  f  f  K{x,t)  bi{x)  bj{t)  dx  dt 
Jo  Jo 

and  {6i,  621  •  • .}  is  a  multi-wavelet  basis  of  £^[0, 1].  This  computation  is  done  for 
the  multi-wavelet  basis  of  order  k  =  A,  for  various  sizes  n. 

In  addition  the  inverse  matrix  (/  —  7’)'"‘  is  obtained  by  the  Schulz  method. 
Table  2.2  displays,  for  various  precisions  c,  the  average  number  of  elements  per 
row  in  the  matrices  I  —  T  and  (I  —  T)~^.  Fig.  2.3  displays  the  matrices  for  n  =  128 
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Figure  2.2:  Functions  are  graphed  for  A:  =  4  (top  graph)  and  k  =  5 

(bottom).  Each  function  (given  in  Table  2.1)  is  a  polynomial  on  the  interval 
(0, 1),  is  an  odd  or  even  function  on  (  —  1,1),  and  is  zero  elsewhere. 


and  €  =  10"^. 
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Table  2.2:  The  average  number  of  elements  per  row  of  the  matrices  I  —  T  and 
(/  —  T)~^,  where  T  is  defined  in  Eq.  (2.26),  is  tabulated  for  various  precisions  e 
and  various  sizes  n.  Here  fc  =  4. 

c  =  10-2  g  ^  10-2  g  _  10-4 


n 

I-T 

1 

1 

I-T 

(I-T)-^ 

I-T 

1 

1 

32 

8.8 

9.7 

19.3 

19.6 

22.8 

23.6 

64 

9.3 

10.0 

25.8 

26.0 

31.9 

32.6 

128 

9.9 

10.1 

29.2 

29.4 

38.2 

38.8 

256 

11.8 

11.8 

30.1 

30.3 

41.9 

42.7 

2.6  Discussion 

The  results  of  the  previous  subsection  demonstrate,  for  a  particular  integral  oper¬ 
ator,  that  the  multi- wavelet  representations  are  sparse.  The  matrix  has  a  peculiar 
structure  in  which  the  non- negligible  elements  are  contained  in  blocks  lying  along 
rays  emanating  from  one  corner  of  the  matrix.  Furthermore,  the  inverse  matrix 
shares  that  structure.  This  property  is  a  general  characteristic  of  integral  opera¬ 
tors  with  non-oscillatory  kernels  that  possess  diagonal  singularities,  and  further 
examples  will  be  given  in  Chapter  4. 

The  kernel  K{x,  t)  =  log  |x— f  |  of  the  previous  subsection  was  chosen,  however, 
because  the  projections  Kij  could  be  computed  analytically,  thereby  avoiding  use 
of  quadratures.  The  difficulty  here  with  quadratures  is  that  they  would  be  re¬ 
quired  for  each  element  Kij  (in  two  dimensions),  and  would  have  to  cope  with  the 
singularity  of  the  logarithm.  It  was  felt  that  the  analytical  computation  would 
be  more  efficient.  In  fact,  the  analytical  computation,  which  requires  integrating 
monomials  (0  <  j  <  k)  against  the  logarithm  and  combining  the  results  with 
large  coefficients,  is  a  very  poorly-conditioned  procedure.  The  computations  de¬ 
scribed  above  required  quadruple-precision  arithmetic  to  obtain  single-precision 
accuracy  for  n  as  small  as  64.  This  procedure  is  not  recommended. 

The  fault  lies,  of  course,  not  with  the  idea  of  projection  to  the  multi-wavelet 
basis,  but  with  the  method  of  projection.  The  integration  should  be  performed 
numerically,  with  quadratures.  As  mentioned  above,  such  a  procedure  would 
require  use  of  quadratures  for  each  matrix  element  Kij,  or  potentially  order 
0(n  log  n)  times.  A  more  efficient  procedure  is  to  use  the  Nystrom  method,  in 
which  only  n  quadrature  applications  are  required.  Numerical  quadratures  and 
a  vector-space  analogue  of  the  multi-wavelet  bases  are  developed  in  Chapters  3 
and  4;  these  tools  enable  efficient  solution  of  integral  equations  using  Nystrom’s 
method. 


2.7.  APPENDIX:  TABULATION  OF  BASIS  FUNCTIONS 


25 


Figure  2.3:  Matrices  re-presenting  the  operators  1  —  JC  (top)  and  (1  —  (bot- 
tom),  with  JC  defined  by  Eq.  (2.25),  expanded  in  the  multi-wavelet  basis  of  order 
k  =  4,  for  n  =  128.  The  dots  represent  elements  above  a  threshold,  which  is 
determined  so  as  to  bound  the  relative  truncation  error  at  e  =  10"^. 

2.7  Appendix:  Tabulation  of  Basis  Functions 

In  the  following  table  we  present  the  functions  /i, . . . ,  /t,  expressed  in  expansions 
of  orthonormal  polynomials  for  the  interval  (0, 1),  and  tabulated  for  k  =  4,S,  12. 
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Table  2.3:  Coefficients  for  the  orthonormal,  vanishing-moment  functions 
fi,...,fk,  for  various  k,  for  argument  x  in  the  interval  (0,1).  Each  function 
is  given  in  the  form  f{x)  =  Cj  •  pj{x),  where  pj{x)  =  Pj{2x  -  1)  y/2j  -f  1, 
the  Legendre  polynomial  of  degree  j,  shifted  to  the  interval  (0,1)  and  normalized. 
The  function  fi  is  extended  to  the  interval  (—1,1)  as  an  odd  or  even  function, 
according  to  the  formula  fi{x)  =  (— l)‘‘''^“^/i(— x)  for  x  G  (  —  1,0),  and  is  zero 
outside  the  interval  (—1,1). 

i  =  4 

/i  fi  fs  A 

Co  .0000000000000000  -.1543033499620919  .0000000000000000  .2156454872944857 

Cl  -.1533929977694741  .2672612419124244  .0878668779193509  -.3735089404169980 

C2  .5940885257860046  -.1725163898355886  -.3403069548648863  .4436221311410142 

C3  .3514675116774037  -.6123724356957945  .6135719910778963  -.3423265984407288 


i  =  8 

/i  /a  fz  A 

Co  .0000000000000000  -.0843312289034136  .0000000000000000  .0742292692888546 

Cl  -.0394514921576236  .1460659731254334  .0258430814883384  -.1285688658170082 

C2  .1527949721095086  -.1616317375321096  -.1000898242190157  .1503045161172336 

C3  -.3013157484854446  .0637484169813985  .2127004256180621  -.1036447629483221 

ca  .2049959665542275  .2168517314659207  -.2558883246467202  -.0669121025632302 

C5  .5288069924640356  -.3995643491174248  -.0103051686103246  .3099911232623849 

C6  .2463744895581894  -.4560908549012228  .4962827560403292  -.2333502563584798 

C7  .0441081091391231  -.1633067225540632  .3635809208901253  -.5381471079157126 

A  A  A  A 

Co  .0000000000000000  -.0871939722326493  .0000000000000000  .1309942892720891 

Cl  -.0213434791213834  .1510243900206985  .0192854509208144  -.2268887645206330 

C2  .0826629391872435  -.1841022341475768  -.0746922302404146  .2899079084461328 

C3  -.1829592311019578  .1663892400581053  .1723717448166373  -.3288050657844448 

c<,  .2692181887254343  -.0520293149261530  -.2992501225820789  .3336898392328423 

C5  -.1968270025380475  -.1836375913703662  .4050432626343038  -.2900996008838325 

C6  -.1805321245713286  .4401344681347232  -.4003612351829521  .1945752095162505 

C7  .5616280477927965  -.4217364916434134  .2245825563246441  -.0763964570434043 
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Table  2.3:  (continued) 


i  =  12 


h 

/a 

Co 

.0000000000000000 

-.0592042739052351 

Cl 

-.0179405992366968 

.1025448104290915 

C2 

.0694836420647080 

-.1233585459417166 

C3 

-.1530096755279723 

.1032398584967429 

C4 

.2201077082230140 

-.0064287469139523 

C5 

-.1431401336550512 

-.1680946588317245 

C6 

-.1805072879652185 

.2749158159062862 

C7 

.3717453015731024 

-.0081064654950250 

ca 

.4413074075372503 

-.4549238393986653 

Co 

.1881229915369518 

-.3636376752880584 

Cio 

.0395553032931457 

-.1281176981334961 

cn 

.0041396058259610 

-.0236611435924639 

/s 

/a 

Co 

.0000000000000000 

-.0490852233087548 

Cl 

-.0104273737468520 

.0850181006716274 

C2 

.0403850448662383 

-.1054116332841276 

C3 

-.0918988661694854 

.1041544568835931 

C4 

.1514809211499543 

-.0626055685259494 

Cs 

-.1754801708549531 

-.0348056650511157 

C6 

.0929109878353303 

.1662387148811842 

C7 

.1216196841160863 

-.2280521421975709 

Ca 

-.2844627881902686 

.0647562432345634 

Cg 

.0147751774362919 

.2752581716809581 

ClO 

.4359353132957034 

-.2219439549151931 

Cll 

.3761307967348913 

-.5043899953764185 

h 

fio 

Co 

.0000000000000000 

-.0663210971285699 

Cl 

-.0090771349437703 

.1148715098403935 

C2 

.0351555924684996 

-.1458392917293345 

C3 

-.0816263402165753 

.1609203666648590 

C4 

.1447831048236192 

-.1500812540366636 

cs 

-.2066628313763642 

.0983488566098737 

C6 

.2287027892690286 

.0060095635684675 

C7 

-.1579200668486585 

-.1566209531327825 

c« 

-.0366328605538905 

.3115846050759926 

Cg 

.295:802756507381 

-.3921104938549786 

Cio 

-.4314237351640850 

.3270789009891432 

Cll 

.2765929192982552 

-.1437437082969879 

/a  f< 

.0000000000000000  .0516869543201819 

.0132868156385064  -.0895244309710467 

-  .05 14596 156920636  . 1096491458548548 

.  1 156254051056246  - .  1016897863958198 

-.1813181292375394  .0408535191494144 

.1765547047788799  .0858956775115070 

-.0058679070010733  -.2211070510482991 

-  .273694 16472 16635  . 185301 1503121983 

.1993616369248909  .1593931614774968 

.4671192165219669  -.3452248228872576 

.2822961661345380  -.4313594000122254 

.0847279654920732  -  .2098 115229292681 

h  /a 

.0000000000000000  .0549972807348283 

.0094614910169722  -.0952580845454926 

-  .0366441971390245  . 1 194424834832239 

.0842044896200951  -.1245947511300136 

-.1439433554135965  .0954606130089946 

.1860088057424346  -.0149295180909320 

-  .1548304378705 184  -.11 54462360082784 

.0001986646990547  .2421338205360575 

.2270020705291812  -.2425421502237302 

-.2763185483542785  .0010478959001006 

-.1129177018240620  .3713168669707026 

.5211 079830677339  -  .4243990923037588 

fll  fl2 

.0000000000000000  .0971353300365081 

.007879838956393 1  - .  1 682433268332032 

-  .0305184850489067  .2165192505829722 

.0716769100690038  -.2529614538466701 

- .  1321728306705 1 1 8  .2777662732296559 

.2074571615484287  -  .2875975640399399 

-.2845256425672342  .2777726278131828 

.3405574275289650  -.2447881538876456 

-  .3475763522283202  . 1 8979446572855 1 1 

.2872585567373639  - .  1218486998229358 

-.1727059110267771  .0578048706132345 

.0572056487761727  -  .0 152658200454800 


Chapter  3 

Numerical  Quadratures 


An  integral  equation  to  be  solved  numerically  must  be  converted  into  a  finite¬ 
dimensional  problem.  The  equation’s  solution,  completely  characterized  by  its 
infinite  expansion  in  a  chosen  basis,  is  approximated  by  a  finite  truncation  of 
that  expansion.  Alternatively,  the  solution  may  be  represented  by  its  values  at  a 
finite  set  of  points. 

In  the  previous  chapter,  a  projection  method  wcis  used,  in  which  the  infinite 
expansion  in  a  multi-wavelet  basis  was  approximated  by  a  truncated  expansion 
containing  a  finite  number  of  terms.  By  contrast,  in  this  chapter  and  the  next,  we 
use  the  method  of  Nystrom,  in  which  the  integral  is  approximated  by  a  quadra¬ 
ture  (a  weighted  average  of  values  of  the  integrand  at  selected  points).  The 
quadrature’s  rate  of  convergence  to  the  integral,  as  the  number  of  points  in¬ 
creases,  affects  the  amount  of  computation  needed  to  achieve  a  given  accuracy 
in  the  solution  to  the  integral  equation.  The  convergence  rate  depends  on  the 
behavior  of  the  integrand  and  in  particular  on  whether  the  integrand  is  singu¬ 
lar  in  the  interval  of  integration.  In  this  chapter  we  develop  rapidly-convergent 
quadratures  for  several  types  of  singularities  encountered  in  physical  problems. 


3.1  Nystrom  Method 

A  linear  Fredholm  integral  equation  of  the  second  kind  may  be  written  in  the 
form 

f{x)-p{x)j  K{x,t)  f{t)dt  =  g{x),  (3.1) 

where  the  function  p  is  the  coefficient,  K  is  the  kernel,  g  is  the  right  hand  side, 
cind  /  is  the  unknown.  The  coefficient  p  is  often  included  as  part  of  the  kernel  A', 
but  we  separate  them,  for  in  the  following  development  we  allow  p  considerable 
freedom  (including  oscillatory  behavior)  while  K  is  more  restricted.  We  use  the 
symbol  /C  to  denote  the  integral  operator  of  Eq.  (3.1),  which  is  given  by  the 
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formula 

(X:/)(x)  =  p{x)  f  K{x,  t)  f{t)  dt, 

Ja 

for  all  /  €  C[a,b]  and  x  €  [a,b].  Then  Eq.  (3.1)  written  in  operator  form  is 

(I-K:)f  =  g.  (3.2) 

The  Nystrom,  or  quadrature,  method  for  numerical  solution  of  integral  equations 
approximates  the  integral  operator  K  by  the  finite-dimensional  operator  R,  char¬ 
acterized  by  weights  wi,W2, . . .  ,Wn  €  Tt  and  points  xi,X2, . .  ■  ,Xn  €  [a,  6],  and 
given  by  the  formula 


{Rf){x)  =  p{x)  ^  Wj  K{x,  Xj)  fixj), 

for  all  /  €  C(a,  5]  and  x  €  [a,  6].  Substitution  of  R  for  K  in  Eq.  (3.2),  plus  the 
requirement  that  the  resulting  equation  hold  for  x  =  xi,X2, . . .  ,Xn,  yields  the 
following  system  of  n  equations  in  the  n  unknowns  /i, /j, . .  • ,  fn' 

n 

fi  -  p{xi)  fi  =  =  1, . . .  ,n. 

i=i 

The  approximation  (/i,...,/n)  to  the  solution  /  of  Eq.  (3.1)  may  be  extended 
to  all  X  €  (a,  b]  by  the  natural  formula 


fnix)  =  gix)  +  p{x)Y,^j  f,,  (3.3) 

»=i 

which  satisfies  //i(x,)  =  /i  for  i  =  1, . . .  ,n.  A  bound  on  the  error  e/j  =  /  —  /r  of 
the  Nystrom  solution  is  given  in  Chapter  2  by  (2.12), 

lM<||(/-;c)-‘ii-||(^- Wr1|. 

The  error  depends,  therefore,  on  the  conditioning  of  the  original  integral  equation, 
as  is  apparent  from  the  term  ||(/  —  ^)~*||,  and  on  the  fidelity  of  the  quadrature 
R  to  the  integral  operator  K.  It  is  not  necessary  that  \\}C  —  7?||  be  small,  rather 
merely  that  R  approximate  K  well  near  the  solution  /.  In  this  chapter  we  develop 
quadrature  rules  that  have  this  property;  nevertheless,  they  are  of  a  somewhat 
different  form  than  R  and  are  defined  only  on  the  mesh  points  xi, . . . ,  x„  rather 
than  the  whole  interval  [a,  6]. 
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3.2  Corrected  Trapezoidal  Rules  for  Singular 
Functions 

The  quadrature  method  we  develop  arises  from  a  method  developed  by  Rokhlin 
[16].  Some  changes  have  been  made  with  the  aim  of  solving  non-periodic  integral 
equations;  in  addition,  one  improvement  permits  higher  order  quadratures  in 
practice.  These  differences  from  [16]  are  noted  as  they  are  presented. 

It  is  well  known  that  the  trapezoidal  rule  for  integration  can  be  modified  at  the 
ends  via  the  Euler-Maclaurin  summation  formula  to  a  rapidly  convergent  rule, 
provided  the  integrand  is  sufficiently  differentiable.  We  will  suppose,  instead, 
that  the  integrand  is  singular  at  one  end  of  the  interval  and  the  form  of  the 
singularity  is  known.  In  this  case  a  modification  at  that  end  may  be  determined 
so  that  the  corrected  trapezoidal  rule  is  rapidly  convergent. 

3.2.1  Differentiable  Integrands 

We  begin  with  the  assumption  that  the  integrand  is  differentiable  throughout  the 
interval  of  integration.  For  positive  integers  /,  n,  and  m  =  2/  -f  2,  and  a  function 
/  6  C’"[a,6],  the  Euler-Maclaurin  summation  formula  is  given  by  the  equation 
(see,  e.g.,  [18]) 

f'  fix)  dx  =  T„(/)  -f-  /?-(/,  a)  -  h)  -H  E'^U).  (3.4) 

Ja 

where 

TM)  =  A  (i/W  + /(<■  +  A)  +  •  •  •  + /(4  -  A)  + 

m/2-1  n 

D7U.X)  =  E  A“(2^/"‘-‘'(x) 

K(f)  =  A"  E  -  0  -  gm 

771. 

ml 

and  h  =  {b  —  a)/n  and  a  <  ^  <  b.  Here  Bi  are  the  Bernoulli  numbers, 

^2  =  -,  =  ^6  =  — ,  Bs  =  -  — 

The  derivatives  which  appear  in  D^{f,a)  and  D'^[f,b)  may  be  approximated 
by  finite  differences,  to  obtain  high-order  quadrature  rules  which  depend  only  on 
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the  values  of  /  at  equispaced  points.  Using  the  Taylor  expansion  for  /  about  the 
point  a  we  have,  for  2  =  1, 2, . . . ,  m  —  2, 


"•-2  (i  k\j 

f{a  +  *  A)  =  f{a)  +  ^ 

i=i  J' 


{i 

(m  —  1)! 


(3.5) 


with  a  <  Vi  <  a  A  ih.  Eqs.  (3.5)  can  be  considered  to  be  a  system  of  equations 
in  the  unknowns  hf^^\a),  and  the  matrix  A”'  of 

this  system,  given  by  the  formula 


A'"  = 


/ 


1 

2 


± 

2! 

2! 

2! 


(m-2)! 


y  m  —  2 


2! 


(m-2)'"-* 
(m-2)'  / 


is  non-singular,  since  the  functions  /j\  for  j  =  1,2, .. .  form  a  Chebyshev  sys¬ 
tem.  We  define  the  vector  u”  =  {u^i,  •  •  • ,  of  finite  differences  by  the 

expression 


/ 


/(a  -K  h)  -  /(a) 
/(a  -b  2h)  -  f{a) 


\ 


and  from  Eqs.  (3.5)  we  obtain,  for  j  =  1, . . . ,  m  -  2,  the  irrcr  bound 


h^f^^\a)\<e^h^-^  sup  |/<”‘-^He)|, 

2)/i] 


(3.6) 


where  tm  is  independent  of  [a,  6],  A,  and  the  function  /. 

The  expressions  for  the  derivatives  are  used  to  define  the  left-end  correc¬ 
tion  iD'^(f,a)  by  the  formula 


m/2-l  D 

L  ^2i  m 

g  *(2.)l"“-- 


(3.7) 


Similarly,  in  Eqs.  (3.5)  we  replace  a  by  6  and  i  by  —i  to  obtain  finite-difference 
expressions  •  •  •  j^^m-2)^  fo*"  ^^e  derivatives  at  the  right  endpoint 

(/i/'(6),  h^/<^)(6), . . . ,  h'"~^/^'""^^(6))^,  and  we  define  the  right-end  correction 
iDjJ*(/, 6)  by  the  formula 


m/2-1 


5. 


,D:(f,b)=  £ 


i=l 


(2i)! 


(3.8) 
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We  thus  obtain  the  corrected  trapezoidal  rule 

iTTC/)  =  T.if)  +  .D:(f,a)  -  rD^if^b)  (3.9) 

that  depends  only  on  the  values  of  /  at  the  equispaced  points  a,  a  +  h, . . . , 
b  —  h,  b.  Combination  of  the  Euler-Maclaurin  formula  (Eq.  3.4)  with  the  finite 
differences  error  bound  (3.6)  yields  the  bound 

.rn/)-  <  ,e„  (4 sup 

where  le^  is  independent  of  [a,  6],  h,  and  the  function  /. 

Alternatively,  the  corrections  to  the  trapezoidal  rule  may  be  concentrated  in 
the  subintervals  [a,  a  +  h)  and  (6  —  h,b]  by  replacing  h  hy  h'  =  h/{m  —  1)  in 
Eqs.  (3.5)  to  define  revised  endpoint  corrections  ■iD'^{f,a)  and  2D^{f,b)  by  the 
formulae 

,D"(/,u)  =  (3.10) 

,D:U,b)  =  (3-11) 

These  endpoint  corrections  give  us  the  “crowded”  corrected  trapezoidal  rule 

2C(/)  =  r„(/)+  2£>:(/,a)-  2D:(/,6).  (3.12) 

For  the  quadrature  2T^  we  obtain  the  error  bound  (as  for  iT^) 

:C(/)-  /‘/(x)<ii  <  sup  |/'”'K)| . 

•'<»  «e[o,6i  '  ' 

where  jCm  is  independent  of  [a,  6],  h,  and  the  function  /.  This  rule  has  the 
advantage  over  the  equispaced  rule  that  the  constant  of  the  error  term  is  smaller; 
it  has  the  disadvantage  that  the  coefficients  are  larger  and  hence  produce  larger 
rcund-off  errors.  Independently  of  these  two  characteristics,  the  crowded  rule 
is  a  suitable  starting  point  for  construction  of  rules  for  integral  equations  with 
singular  kernels,  as  we  shall  see  below. 

Coefficients  of  finite-difference  expressions  for  the  derivatives,  equispaced  cor¬ 
rected  trapezoidal  rules,  and  crowded  corrected  trapezoidal  rules  are  given  in  the 
chapter  appendix. 
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3.2.2  Singular  Integrands 

We  now  consider  integrands  of  the  form 

f{x)  =  0(ar)  six)  +  ^(i)  (3.13) 

for  all  X  €  (0,6],  where  <^,0  6  C"‘[0,6]  and  the  function  s  contains  the  singular 
part  of  /.  In  particular,  we  construct  quadratures  for  /,  where  s  is  given  by  the 
formula  s(i)  =  log(x)  or  the  formula  s(x)  =  x“  with  0  <  |a|  <  1. 

We  define  to  be  the  trapezoidal  rule  on  [0, 6]  minus  the  left-end,  by  the 
formula 

Ku,  4)  =  /.  (/(A)  +  ••+/(('-  M  +  Im) . 

where  h  =  b/n,  and  define  the  right-end  corrected  rule  T"  by  the  formula 

=  2D::if,b),  (3.14) 

where  the  term  6)  is  the  finite-difference  approximation  to  the  Euler- 

Maclaurin  correction,  as  defined  in  Eq.  (3.11).  Given  a  positive  integer  /  and  two 
finite  sequences  =  (/3t,  /?2,  •  •  • ,  A)  and  x  =  (Xi,  X2,  •  ■  • ,  X/),  with  0  <  Xi  <  X2  < 
•  •  •  <  X/  ^  we  define  the  left-end  correction  L^  by  the  formula 

1=1 

We  now  define  the  linear  mapping  :  £'(0,  6]  — >  7^  by  the  formula 

rf(/,4)  =  r"(/,6)  +  i';(/).  (3.15) 

Although  L^  and  depend  on  x>  this  sequence  is  generally  implied  (see  Remark 
3.1  below)  and  we  omit  the  symbol  x  from  our  notation.  The  mapping  with 
suitable  choice  of  and  x  will  be  used  ais  a  quadrature  formula  for  functions  / 
of  the  form  g'wen  in  Eq.  (3.13).  We  now  show  how  /3  may  be  chosen,  given  x- 
For  positive  integers  k  and  n  we  consider  the  following  system  of  linear  equa¬ 
tions  with  respect  to  the  unknowns  =  (/?", . . . ,  02k)' 

E0J{X,hY  =  /‘x'</x-r:(x',6) 

Jo 

(3.16) 

E/JJiXihf^X.h)  =  f'’x’sU)dx-T:{x's{x),b), 

J  =  l 

for  i  =  0, 1, . . . ,  A:  —  1.  We  will  see  below  that  this  system  has  a  unique  solution 
d"  and  for  functions  /  of  the  form  given  in  Eq.  (3.13)  the  rule  T^”  is  a  quadrature 
with  convergence  of  order  at  least  A:  —  1. 
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Remark  3.1  In  Rokhlin  [16],  the  points  xii  •  •  •  ?  X2ik  are  chosen  to  be  equispaced, 
namely  x,  =  if  {2k).  While  this  choice  leads  to  satisfactory  quadratures  for  k  <3, 
larger  values  of  k  correspond  to  coefficients  ,  02k  which  are  large  enough  to 

introduce  substantial  round-off  errors  to  fixed-precision  computations.  We  have 
found  that  alternative  spacing  of  the  points  can  be  chosen  so  as  to  delay  the 
growth  in  the  coefficients.  Letting  Xi»  •  •  • » X2fc  be  half-Chebyshev  points 

results  in  . . . ,  02k  which  are  satisfactory  for  k  <  5.  Coefficients  for  both 
equispaced  points  and  half-Chebyshev  points  are  given  in  the  chapter  appendix. 

Our  purpose  in  this  chapter  is  the  development  of  quadratures  for  integral 
operators  with  diagonally-singular  kernels.  We  consider  kernels  of  the  form 

K(x,  t)  =  <^(i,  t)  s(|x  —  t\)  -f  0(z,  f), 

where  €  C"*([a,6]  x  [a,  6])  and  s  is  as  introduced  above.  To  compute  the 
value 

{Kf){x)  =  pix)  f'‘K{x,t)mdt,  (3.18) 

Ja 

the  interval  [a,  6]  is  divided  into  [a,i]  and  [i,  6).  The  integrand  is  singular  on  one 
end  of  each  of  these  intervals,  so  the  quadrature  is  applicable.  We  augment 
our  notation  by  defining  the  quadrature  cTn"  for  an  interval  [c,  6]  with  a  left-end 
singularity  by  the  formula 

trr(/(i),c,(.)  =  rf  (/(x  +  c),6-  c). 

Similarly,  we  define  the  quadrature  rT^”  for  an  interval  [a,  c]  with  a  right-end 
singularity  by  the  formula 

Rrr(/(x),a,c)  =  rr(/(c-x),c-a). 

Now  we  define  a  quadrature  for  the  integral  in  Eq.  (3.18)  by  defining  the  mapping 
•  ^*([o,  b]  X  [a,  6])  — ♦  by  the  formula 

r‘(G)  =  nTr{Gi,a,xi)  +  irt7(G„Xi,6),  (3.19) 

and  further  define  r^„(G)  =  cTr{Go,<i.b)  and  r‘„(G)  =  (G„,<i,6).  Here 

G  :  [a,  6]  x  [a,f>]  — »  TZ  and  G,  :  [a,  6]  — ♦  7?.  is  the  restriction  of  G  defined  by  the 
formula  Gi{t)  =  G{xi,  t)  where  Xi  =  a  +  i  h  for  t  =  0, 1, . . . ,  n,  and  h  =  {b  —  a)/n. 
The  quadrature  Tj^-  has  convergence  of  order  k  —  I,  uniformly  in  z,  as  we  will 
show  in  the  next  section. 
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3.3  Analytical  Properties  of  the  Corrected  Rules 

The  following  lemma  is  a  restatement  of  a  classical  result  (see,  e.g.,  [11]). 

Lemma  3.2  Suppose  that  the  function  s  :  (0, 6]  -+  72.  is  given  by  the  formula 
s(x)  =  log(x)  or  the  formula  s(x)  =  x"  with  0  <  [a]  <  1.  Then  the  system  of 
Eqs.  (3.16)  has  a  unique  solution  . . . , 

The  next  two  lemmas  are  proven  in  [16].  Lemma  3.3  states  that  the  coefficients 
have  limiting  values  as  n  oo  and  that  the  rate  of  convergence  to 
these  values  depends  on  the  difference  between  the  orders  of  correction  k  on  the 
left  end  and  m  on  the  right  end.  Lemmas  3.3,  3.4,  and  3.5  will  be  used  in  the 
proof  of  Theorem  3.6,  below. 

Lemma  3.3  (Rokhlin)  Suppose  that  k  and  m  are  positive  integers  with  k  <  m, 
the  function  s  is  as  specified  in  Lemma  3.2,  and  for  each  positive  integer  n,  the 
coefficients  =  (^[‘, . . . ,  solution  to  the  linear  system  of  Eqs.  (3.16). 

Then  there  exist  coefficients  fi  =  {fii,  ■  ■  • ,  fiik)  end  a  constant  c  >  0  such  that 

for  i  =  I,. . .  ,2k  and  all  n. 


Lemma  3.4  (Rokhlin)  Suppose  that  the  function  s  is  one  of  the  singular  func¬ 
tions  specified  in  Lemma  3.2,  k  is  a  positive  integer,  and  ip  €  C*[0, 6]  satisfies 

0  =  (,3(0)  =  <p\0)  =  .  • .  =  (,?‘*'Ho)- 

Then  the  function  w  =  ip  •  s  is  defined  on  the  closed  interval  [0, 6]  and 

0  =  u;(0)  =  «;'(0)  =  •  •  •  =  u;^*'”*^(0). 


For  the  function  w  —  ip  ■  s,  whose  A:th  derivative  is  unbounded  near  0,  the 
simplified  Euler-Maclaurin  error  expression 

with  0  <  ^  <  6,  is  not  useful.  In  this  case  we  substitute  the  following  bound. 

Lemma  3.5  Suppose  that  k  and  n  are  positive  integers  and  the  function  w  : 
[0,6]  — >  72  is  as  specified  in  Lemma  3.4,  above.  Then  the  error  El^{w)  of  the 
Euler-Maclaurin  formula  (Eq.  3.4)  has  the  bound 


E'^{w)  <  h^  sup 

x6[0,l] 


Bkjx)  —  Bk 
k\ 


(3.20) 
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Thus,  there  is  a  constant  Ck  >  0,  independent  of  n,  b,  and  the  function  w,  such 
that 

rf{xv,b)-  j\{x)dx^^  <  hJ^Ck  J^\w^'‘\x)\dx,  (3.21) 

where  h  =  b/n. 


Proof.  The  bound  (3.20)  is  immediate  from  the  definition  of  E!^  given  in 
Eq.  (3.4).  The  error  bound  (3.21)  then  follows  from  the  combination  of  the 
Euler-Maclaurin  formula,  the  observation  that  D^{w,0)  =  0,  the  definition  of  T'k 
(Eq.  3.14),  and  bound  (3.6).  D 

The  following  theorem  is  the  foundation  for  the  corrected  trapezoidal  rules 
for  singular  functions.  It  is  a  slight  generalization  of  a  theorem  found  in  Rokhlin 
[16],  in  that  the  error  bound  is  established  for  integrals  taken  on  subintervals 
[0,  ^6],  for  1  =  1,...  ,  n,  of  the  interval  [0, 6],  provided  that  the  trapezoidal  points’ 
spacing  h  =  b/n  is  preserved. 


Theorem  3.6  Suppose  that  k  and  m  are  positive  integers  with  2  <  k  <  m  and 
the  function  f  :  (0, 6]  — »  is  given  by  the  formula  f  =  <f>  •  s  +  ip,  where 
<f>,ip  €  C”*[0, 6]  and  the  function  s  is  as  specified  in  Lemma  3.2.  Further  suppose 
that  for  each  positive  integer  n,  is  the  solution  of  Eqs.  (3.16)  and  T^"  is 
defined  by  Eq.  (3.15).  Then  there  exists  c  >  0  independent  of  b  and  the  function 
f  such  that 


rf  (/,  i/n  ■  b)  -  If"  /(i)  dx  <cb  ft‘-‘  (3-22) 

for  all  n  and  i  =  1, . . . ,  n. 


«€(o,61 


Proof.  vVe  write  the  function  /  as  a  sum  of  two  functions;  one  is  integrated 
exactly  by  the  quadrature  and  the  other  has  several  zero  derivatives  at  0. 

Let  P(v>)  denote  the  fc-term  Taylor  expansion  of  a  function  €  C’"[0,  b]  about 

PMM  =  E  x'. 

j=0  J- 

Now  we  define  pr  =  d>  —  P{<f>)  and  ip^  =  ip  —  P{ip)  so  we  have 

0  =  .^,(0)  =  0p(O)  =  . . .  =  «i<*-*>(0)  =  T^('’-*)(0).  (3.23) 


We  further  define  fp  =  P{p)  •  s  +  P{ip)  and  /^  =  •  s  +  ipr,  so  that  f  =  fp  + 

and  we  let  bi  =  j^b.  Now  we  bound  the  error  by  the  inequality 


TfUM-t!(x)dx 

Jo 


/‘'/,(x)<ix 
Jo 

Tf'i/rM-  t  Ux)dx 
Jo 
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where  by  Eqs.  (3.16)  the  first  term  on  the  right  vanishes.  By  the  definition  of 
jf  (Eq.  3.15),  the  second  term  satisfies 

Tf'Ur.b,)-  t  fr{x)dx  <  T;y,.bi)-  I'"  Mx)dx  ■ 

JO  JO 

By  Eq.  (3.‘23)  and  Lemma  3.4  the  function  fr  =  <i>T’S  +  '4>T  vanishes  at  0  and  h«is 
vanishing  derivatives, 

o=/,(o)  =  /:(o)  =  ...  =  /f-«(o), 

from  which,  in  combination  with  error  bound  (3.21),  we  obtain 

TlVr,  bi)  -  Mx)  dx  <  Cfc-  |/(*'-*)(x)|  dx  (3.24) 

for  some  constant  Cfc  independent  of  n.  6,  i,  and  the  function  /,.  For  the  remaining 
term,  we  define  =  A;!  sup  |<^^*')(x)|  and  M,/,  =  Arlsup  both  suprema 

taken  for  x  €  [0, 6],  therefore  finite,  and  =  sup  taken  over  positive  integers 
i  and  j  =  1, . . . ,  2A:.  Mjs  is  finite  by  Lemma  3.3,  and  we  obtain 

<  M«£(M<,(x,A)*WXiA)l  +  M»(xiA)*) 

J=l 

<  M0  {2k)  {M^  +  M^)  h'^-\  (3.25) 

for  sufficiently  large  n.  Combining  bounds  (3.24)  and  (3.25)  yields  (3.22).  □ 
Note  that  Theorem  3.6  establishes  that  the  quadrature  T^''{f,b)  converges 
to  the  integral  of  /  on  the  fixed  interval  [0,6],  with  order  of  convergence  at  least 
A:  —  1.  Additionally,  however,  it  establishes  the  same  convergence  on  subintervals 
[0,  ^6]  of  [0, 6]  with  correspondingly  fewer  quadrature  points.  This  characteristic 
is  essential  for  proper  treatment  of  non-periodic  integral  operators;  in  fact,  it 
assures  uniform  convergence  of  the  quadratures  T^-  defined  in  Eq.  (3.19). 

Corollary  3.7  Suppose  that  k  and  m  are  positive  integers  with  k  <  m  and  that 
the  kernel  K  :  [a,  6]  x  [a,  6]  — +  77  is  given  by  the  formula 

K{x,t)  =  <l>{x,t)s{\x  -  /|)  +  0(x,Z) 

where  V’  €  C’"([a,6]  x  [a,  6])  and  s  is  as  specified  in  Lemma  3.2.  We  further 
suppose  f  €  ^"[a,  6]  and  define  the  function  G  :  [a,  6]  x  [a,  b]  TZ  by  the  formula 
G{x,t)  =  K{x,t)  f{t).  There  exists  c*  >  0  such  that 
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Proof.  We  observe  that  and  defined  by  the  formulae 


,v,p 

x.t6[«i.6J  dt’‘ 


are  finite,  then  apply  Theorem  3.6.  □ 

The  quadratures  T^^  achieve  uniform  convergence  and  thus,  combined  with 
the  kernel  K  and  the  coefficient  p,  represent  an  operator  R  which  approximates 
the  integral  operator  K.  defined  above.  An  issue  that  arises  in  computing  Rf 
for  a  function  /,  however,  is  that  T^-  requires  the  values  of  /  at  non-equispaced 
points  in  the  interval  [a,  6]  in  addition  to  the  points  a,a  +  h, . . .  ,b  —  h,b.  This 
issue  can  be  handled  by  using  A:-point  interpolation,  for  /  6  C"‘[a,  6].  We -have 
the  estimate 


no -t.  <=<(()  n<‘+{i+k)i>) 


t=i 


where  ci(^), . . .  ,Cfc(^)  are  Lagrange  interpolation  coefficients  and  is  chosen 
such  that  ^  €  [a  +  (1  +  j^)  h,a  +  {k  +  j^)  h].  Replacing  f{()  in  T^iiNf)  by  the 
interpolation  preserves  our  error  estimates. 


3.4  Numerical  Examples 

In  this  section  we  present  numerical  examples  of  the  corrected  quadrature  rules 
applied  to  differentiable  and  singular  integrands.  FORTRAN  routines  were  writ¬ 
ten  which  incorporate  the  corrections  developed  in  Section  3.2,  and  the  quadra¬ 
tures  were  computed  in  double-precision  arithmetic  on  a  Sun  Sparcstation  1.  The 
correction  weights  themselves,  which  are  given  in  the  chapter  appendix,  were  com¬ 
puted  exactly  with  Maple,  and  in  quadruple-precision  arithmetic  in  FORTRAN 
running  on  a  DEC  micro  VAX. 

Table  3.1  shows  the  errors  in  using  quadratures  of  order  4,  8,  and  12  to  ap¬ 
proximate  the  integral  of  a  smooth  function,  comparing  the  equispaced  corrected 
trapezoidal  rule  to  the  “crowded”  corrected  trapezoidal  rule.  Errors  from  the 
uncorrected  trapezoidal  rule  are  also  shown  for  comparison.  We  make  several 
observations: 

1.  Full  single  or  double  precision  accuracy  is  easily  achievable  with  the  higher- 
order  rules. 

2.  The  observed  rate  of  convergence  matches  that  expected  quite  closely,  as 
can  be  seen  by  comparing  the  errors  for  various  n. 
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Table  3.1:  The  quadrature  rules  iT^  of  Eq.  (3.9)  and  2T(D  of  Eq.  (3.12)  are  used 
to  compute  /o[cos(21x)  +  sin(22x)]<ix  and  relative  errors  are  shown  for  various 
values  ofn,m. 


Trapez.  Equispaced  rule  iT^ 


n 

Rule 

rrt  =  4 

m  =  8 

m  =  12 

10 

-0.397E-I-00 

-0.235E-I-00 

0.128E+01 

-0.499E-I-01 

20 

-0.936E-01 

0.840E-02 

-0.119E-01 

0.385E-02 

40 

-0.231E-01 

0.137E-02 

-0.181E-04 

-0.226E-05 

80 

-0.575E-02 

0.108E-03 

0.863E-07 

-0.847E-10 

160 

-0.144E-02 

0.733E-05 

0.617E-09 

0.676E-13 

320 

-0.359E-03 

0.474E-06 

0.286E-11 

-0.189E-15 

640 

•-0.897E-04 

0.301E-07 

0.966E-14 

-0.227E-14 

1280 

-0.224E-04 

0.190E-08 

-0.170E-14 

-0.151E-14 

“Crowded”  rule  2 

'pm 
^  n 

n 

m  =  4 

m  =  8 

m  =  12 

10  ■ 

0.425E-02 

0.103E-04 

0.119E-06 

20 

0.133E-02 

0.866E-06 

0.653E-09 

40 

0.109E-03 

0.462E-08 

-0.420E-11 

80 

0.748E-05 

0.201E-10 

0.121E-11 

160 

0.486E-06 

0.793E-13 

0.810E-12 

320 

0.309E-07 

O.OOOE-l-00 

0.863E-13 

640 

0.195E-08 

-0.227E-14 

-0.172E-12 

1280 

0.122E-09 

-0.170E-14 

0.130E-12 

3.  The  “crowded”  rules  get  a  substantial  jump  on  the  equispaced  rules  and 
show  small  errors  almost  as  soon  as  the  integrand  is  resolved  by  the  quadra¬ 
ture  points. 

4.  The  crowded  rules  ultimately  achieve  somewhat  less  precision  than  the 
equispaced  rules,  due  to  roundoff  error  resulting  from  the  large  correction 
weights,  but  this  effect  is  important  only  if  full  double-precision  accuracy 
is  required. 

In  summary,  the  corrected  rules  should  be  preferred  to  the  trapezoidal  rule  when¬ 
ever  high  accuracy  is  desirable.  The  “crowded”  corrected  trapezoidal  rule  of  order 
12  performs  very  well. 

Table  3.2  shows  the  errors  from  using  the  corrected  trapezoidal  rule  for  singu¬ 
lar  functions  to  compute  the  integral  of  an  integrand  with  logarithm  singularity. 
The  rule  Wtis  applied  for  various  numbers  of  quadrature  points  n  and  various  k, 
and  for  equispaced  and  half-Chebyshev  Xi>  •  •  ••>X2k-  In  each  case  the  correction 
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Table  3.2;  The  quadrature  rule  T^  of  Eq.  (3.15)  used  to  compute  /o  [cos(21x)  + 
sin(22i)+log(i)(cos(23x)+sin(24x))]  dx  and  relative  errors  are  shown  for  various 
values  of  n,k.  The  right  end  is  corrected  to  eighth  order  (m  =  8). 

Equispaced  points  Xi » •  •  •  >  X2k 


n 

A?  =  2 

k  =  3 

II 

k  =  5 

10  ■ 

-0.548E-I-00 

-0.219E-01 

0.429E-01 

-0.608E-02 

20 

0.386E-02 

-0.657E-02 

0.506E-03 

0.394E-04 

40 

0.403E-02 

-0.242E-03 

-0.169E-05 

0.614E-06 

80 

0.663E-03 

-0.341E-05 

-0.236E-06 

0.482E-08 

160 

0.964E-04 

0.381E-06 

-0.971E-08 

0.809E-11 

320 

0.139E-04 

0.545E-07 

-0.363E-09 

-0.593E-12 

640 

0.202E-05 

0.509E-08 

-0.135E-10 

0.258E-14 

1280 

0.292E-06 

0.414E-09 

-0.497E-12 

0.861E-15 

Half-Chebyshev 

points  Xi,--.,X2Jt 

n 

k  =  2 

Jfc  =  3 

k  =  4 

k  =  5 

10 

-0.510E-01 

-0.667E-02 

-0.503E-02 

-O.llOE-02 

20 

0.143E-02 

0.599E-04 

-0.166E-04 

0.553E-05 

40 

0.484E-03 

0.615E-05 

0.732E~06 

0.877E-07 

80 

0.750E-04 

0.349E-06 

0.351E-07 

0.653E-09 

160 

0.108E-04 

0.203E-07 

0.132E-08 

-0.312E-12 

320 

0.156E-05 

0.125E-08 

0.492E-10 

-0.148E-12 

640 

0.227E-06 

0.796E-10 

0.183E-11 

-0.445E-14 

1280 

0.327E-07 

0.519E-11 

0.670E-13 

-O.lOOE-14 

weights  j3  =  (y3i, . . ,  ,^2fc)  were  chosen  to  be  the  limiting  values  (of  Lemma  3.3), 
given  in  the  chapter  appendix.  Note  that  Theorem  3,6  tissures  us  of  convergence 
of  order  at  least  k  —  We  observe: 

1.  Although  we  have  proven  that  the  order  of  convergence  is  at  least  k  —  I, 
the  actual  order  of  convergence  appears  to  lie  between  k  and  +  1  and  the 
convergence  pattern  is  somewhat  irregular. 

2.  The  correction  weights  from  equispaced  points  i  X2k  are  much  larger 

than  those  from  half-Chebyshev  points;  nevertheless,  both  versions  of  the 
quadrature  rule  perform  well  and  lead  to  nearly  full  precision  accuracy  in 
practice. 

The  uncorrected  trapezoidal  rule,  with  the  left  end  omitted,  gives  very  slow 
convergence  for  this  problem  and  is  not  practical  for  achieving  high  accuracy. 

Our  final  examples  of  the  corrected  trapezoidal  rules,  in  Table  3.3,  demon¬ 
strate  the  uniform  convergence  of  the  quadratures  for  various  integral  operators. 
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Table  3.3:  ihe  quadrature  rule  T^^  of  Eq.  (3.19)  is  used  to  compute  F{x)  = 
/(J[cos(21xi)  +  sin(22xf)  +  s(|x  —  f|)(cos(23xf)  +  s\n{2Axt))]  dt ,  with  x  =  ifn  for 
i  =  0, 1, . . . ,  n.  The  function  s  is  singular  and  has  one  of  three  forms  given  below. 
The  relative  C?  errors  for  the  quadratures  with  k  =  A  are  shown. 


Three  Choices  of  Integrand  /(x,  t) ' 


n 

s(x)  =  log(x) 

s(x)  = 

s(x)  =  x'/^ 

10  ■ 

0.302E-03 

0.134E-03 

0.482E-04 

20 

0.527E-05 

0.187E-05 

0.579E-06 

40 

0.681E-07 

0.487E-07 

0.206E-07 

80 

0.573E-08 

0.266E-08 

0.611E-09 

160 

0.283E-09 

0.132E-09 

0.181E-10 

320 

0.118E-10 

0.613E-11 

0.547E-12 

640 

0.454E-12 

0.300E-12 

0.169E-13 

1280 

0.216E-13 

0.155E-12 

0.271E-14 

In  these  examples,  we  approximate  the  integral  Fi  =  /q  f(xi,t)  dt  for  Xi  =  ifn 
and  i  —  0,1,..,, n,  by  the  quadrature  T*,{/),  defined  in  Eq.  (3.19),  The  table 
shows  the  relative  error  of  the  approximations,  defined  by  the  formula 


\  ) 


1/2 


(3.27) 


By  taking  ratios  of  the  errors  for  different  n,  it  can  be  seen  that  (for  fc  =  4)  the 
rate  of  convergence  is  uniformly  of  order  at  leeist  4.  For  the  logarithm  singularity, 
the  order  of  convergence  is  about  4.5,  and  for  the  square  root  singularity,  it  is 
nearly  5.  In  each  case,  n  =  20  produces  roughly  single-precision  accuracy. 


3.5  Appendix:  Quadrature  Weights 

This  appendix  contains  various  weights  and  correction  coefficients  used  in  the 
quadratures  of  the  chapter.  The  coefficients  in  the  finite  difference  expressions  for 
odd-numbered  derivatives,  which  appear  in  the  Euler-Maclaurin  formula,  as  well 
as  the  corresponding  correction  coefficients  are  tabulated  (Tables  3.4  and  3.5  and 
3.6).  The  limiting  values  of  the  correction  weights  used  for  the  singularities  s(x)  = 
log(x)  and  s(x)  =  x"  for  a  =  are  shown  for  equispaced  correction  points 
(Table  3.7)  and  for  half-Chebyshev  correction  points.  The  correction  weights  for 
finite  n,  which  are  very  numerous,  are  not  shown. 
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Table  3.4:  Finite-difference  coefficients  for  odd-numbered  derivatives  to  various 
orders  of  approximation  are  shown.  The  table  entries  are  coefficients  from  the 
formula  =  {l/d)J2^^Cif{x  +  ih)  +  for  j  =  1,3, 5, 7,9  and 

m  =  4,6,8, 10, 12. 


hf^^\x)  h^f^^Hx) 

0{h*)  0{h^)  O(h^)  Q(ftiO)  O(h^)  O(/t^0)  0{h^^) 


d 

2 

12 

60 

840  2520 

2 

6 

288 

Co 

-3 

-25 

-147  - 

-2283  -7381 

-7 

-81 

-8591 

Cl 

4 

48 

360 

6720  25200 

40 

575 

72492 

C2 

-1 

-36 

-450  - 

11760  -56700 

-95 

-1790 

-278313 

C3 

16 

400 

15680  100800 

120 

3195 

640752 

C4 

-3 

-225  - 

14700  -132300 

-85 

-3580 

-979878 

cs 

72 

9408  127008 

32 

2581 

1039656 

Ce 

-10  - 

-3920  -88200 

-5 

-1170 

-774402 

Cj 

960  43200 

305 

399408 

Cs 

-105  -14175 

-35 

-136347 

C9 

2800 

27788 

C\0 

-252 

-2565 

0(/i«) 

h^f^^Hx) 
0(/i«)  0(/iiO) 

0(/i‘2) 

0{h^°)  0(/ll2) 

d 

2 

8 

240 

30240 

2 

24 

2 

Co 

-5 

-49 

-2403 

-420475 

-9 

-605 

-11 

Cl 

18 

232 

13960 

2876868 

70 

5628 

108 

C2 

-24 

-461 

-36706 

-9389763 

-238  ■ 

-23583 

-477 

C3 

14 

496 

57384 

19227792 

462 

58632 

1248 

C4 

-3 

-307 

-58280 

-27098442 

-560  - 

-95802 

-2142 

Cs 

104 

39128 

27147960 

434 

107520 

2520 

Cs 

-15 

-16830 

-19395138 

-210  - 

-83958 

-2058 

Cj 

4216 

9693648 

58 

45048 

1152 

Cs 

-469 

-3229227 

-7  - 

-15897 

-423 

Cg 

645412 

3332 

92 

Cio 

-58635 

-315 

-9 
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Table  3.5:  Equispaced  endpoint  corrections  transform  the  trapezoidal  rule 
into  a  high-order  quadrature  for  functions  with  several  continuous  deriva¬ 
tives.  The  quadrature  rules  are  given  by  the  formula  f{x)  dx  =  Tn{f)  + 
ih/d)ZT=o^  Ci[f{a  +  ih)  +  f{b-  ih)]  +  where  h  =  {b  -  a)/n  and 

TnU)  =  h  [\f{a)  +  /(a  +  A)  +  . . .  +  /(6  -  A)  +  \f{b)]. 


Ojh^)  Ojh*)  Ojh^)  O(h^)  0(/t»°) _ 0(^ 


d 

1 

24 

1440 

120960 

7257600 

958003200 

Co 

0 

-3 

-245 

-23681 

-1546047 

-216254335 

Cl 

4 

462 

55688 

4274870 

679543284 

C2 

-1 

-336 

-66109 

-6996434 

-1412947389 

C3 

146 

57024 

9005886 

2415881496 

C4 

-27 

-31523 

-8277760 

-3103579086 

Cs 

9976 

5232322 

2939942400 

C6 

-1375 

-2161710 

-2023224114 

C7 

526154 

984515304 

Cs 

-57281 

-321455811 

Cg 

63253516 

ClO 

-5675265 
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Table  3.6:  Corrections  crowded  into  the  intervals  [a,  a  +  /i)  and  {b  —  h,  6]  con¬ 
vert  the  trapezoidal  rule  into  a  high-order  quadrature  for  functions  with  sev¬ 
eral  continuous  derivatives.  The  quadrature  rules  have  the  form  /(x)  dx  = 
T4f)  +  Wd)  Er=3'  ^  !/(<■  +  ^A)  +  /(»-  “Aere  h  =  {b-  a)/n 

and  Tn{f)  =  h  [5/(3)  +  /(a  +  A)H - ^  f{b—h)  +  ^f{b)].  These  rules  have  smaller 

error  constants,  but  larger  coefficients,  than  the  equispaced  corrected  trapezoidal 
rules. 


0(A=) 

0(A«) 

0(A») 

0(h^°) 

d 

1 

a 

288 

17280 

89600 

87091200 

Co 

0 

-3 

-125 

-7889 

-41943 

-41374135 

Cl 

4 

30 

13832 

-372570 

4717178004 

C2 

-1 

240 

-57421 

2898654 

-43825028709 

C3 

-190 

133056 

-9112466 

180245487576 

C4 

45 

-130067 

15663360 

-428859839166 

cs 

58744 

-15657582 

649926182400 

C6 

-10255 

9131810 

-649910688834 

C7 

-2897574 

428850243624 

cs 

388311 

-180220260891 

Cg 

43821791596 

Cio 

-4703691465 
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Table  3.7:  The  limiting  values  of  the  correction  weights,  at  equispaced  points,  for 
the  corrected  trapezoidal  rules  for  singular  functions  are  tabulated.  The  corrected 
quadrature  rules  are  given  by  the  formula  f{x)dx  =  T^Xf)  +  ^!Ci=i  + 

+  (h/d)  a  fib  -  h),  where  h  =  (b  -  a)/n,  T^if)  =  h  [/(a  +  /i)  + 

f{a  +  2h)  +  •  •  •  +  fib  —  h)  +  \fib)],  and  the  coefficients  d  and  cq,...,  Cm-2  are 
given  by  Table  3.6. 


01 

02 

03 

04 

0S 

06 


Singularity  5(1)  =  log(i  —  a) 
k  =  2  it  =  3 


0.1601298415357170E+01 

-0.3382558521919485E+01 

0.3627888464434128E+01 

-0.1346628357871812E+01 


0.2228766018460087E+01 

-0.1231212070062612E+02 

0.3157965997308673E+02 

-0.3840391590010434E+02 

0.2267350459115253E+02 

-0.5265893981968890E+01 


k  =  4 

01  0.3093483401777122E+01 

02  -0.3101788376740780E+02 

03  0.13620S9155903270E+03 

04  -0.3 1474 74808. 242 13E+03 

0i  0.4215054 127612634E+03 

06  -0.3287854038787327E+03 

07  0.138801 1671370668E+03 

06  -0.2455521037187227E+02 

09 
010 


k  =  5 

0.4335882754006495E+01 
-0.6931374375731333E-I-02 
0.4590379369965t>28E+03 
-0.1630223863113166E+94 
0.3509093316535128E+04 
-0.4819809483739626E+04 
0.4271587513841620E+04 
-0.2374001 155274844E+04 
0.7547424 129498612E+03 
-0. 1049488 171922296E+03 


Singularity  s(x)  =  {x  -  a) 
k.=  2  k  =  3 

01  0.3338954623777353E+01  0.51568623842001 15E+01 

02  -0.1036918555513964E+02  -0.3718025418572306E+02 

03  0.1238817390561390E+02  0.1050278047880726E+03 

04  -0.4857942974251604E+01  -0.1388205130918806E+03 

0s  0.8799431664143500E+02 

06  -0.2167821653610409E+02 
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Table  3.7:  (continued) 

- 

Singularity  3(x)  =  (x 

—  (continued) 

k  =  4 

k  =  5 

01 

0.7889576157976986E+01 

0.1191640221424121E+02 

02 

-0.1014839102693306E+03 

-0.243121 1398143441E+03 

03 

0.4982052353339497E+03 

0.1823927314310889E+04 

0A 

-0.1241 77860454341  lE+04 

-0.7083666085060807E+04 

05 

0.1751093993580452E+04 

0.1635686640624526E+05 

06 

-0.1419085152097947E+04 

-0.2381679753707451E+05 

07 

0.6179863268019096E+03 

0.2217597560756524E+05 

0S 

-0.1 123274649636003E+03 

-0.1285030321889555E+05 

09 

0.4231931207608981E+04 

010 

-0.6062289570994023E+03 

Singularity  s(z)  =  (i  - 

k  =  2 

Jb  =  3 

01 

0.1076226369733505E+01 

0.140373389574362lE-f-01 

02 

-0.1472473730142097E+01 

-0.6106269754659710E+01 

03 

0.1382935017750346E+01 

0.145821441 1894674E+02 

04 

-0.4866876573417537E+00 

-0.1639615278631818E+02 

05 

0.8952282755716858E+01 

06 

-0.1935738229429337E+01 

k  =  4 

k  =  5 

01 

0.1761334695584808E+01 

0.2301 960768321 394E+01 

02 

-0.1382118344852977E+02 

-0.285624475262 1 036E+02 

03 

0.5459150117813370E+02 

0.1642204330046712E+03 

04 

-0.1173574845498706E+03 

-0.5199121092937581E+03 

05 

0.1507790199321616E+03 

0.101 1593120489479E+04 

06 

-0. 1 14778491 1579322E+03 

-0.1264688737067177E+04 

07 

0.4762309598361213E+02 

0.1027260908792144E+04 

06 

-0.8297842633159577E+01 

-0.5280676861522366E+03 

09 

0.1570202164406574E+03 

010 

-0.2066565945589198E+02 

• 
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Table  3.8:  The  limiting  values  of  the  correction  weights,  at  half-Chebyshev 
points,  for  the  corrected  trapezoidal  rules  for  singular  functions  are  tabulated. 
The  quadratures  are  given  by  faf{x)dx  =  T^(f)  +  /3i  f(a  +  + 

(*/<^)  Ci  f{b  -  h),  where  h  =  (b  -  a)/n,  Xi  =  I  -  cos((2z  -  l)7r/(8fc)), 

T^(/)  =  h  [f{a  A  h)  +  f{a  +  2h)  H - f-  f{b  —  h)  +  lf{b)],  and  the  coefficients  d 

and  Co,. . .  ,Cm-2  are  given  by  Table  3.6. 


/3i 

02 

03 

0A 

05 

0e 


Singularity  s(i)  =  log(i  -  a) 
k=2  k=3 

0.5603992216960789E-01  0.2395413023703419E-01 

0.2870630583396921E+00  0.1743092878178258E+00 

0.2563533589170907E+00  -0.1631724711422584E+00 

-0.9945633942639063E-01  0.9387069242218954E+00 

-0.5880276430135620E+00 


0.1142297718790650E+00 


01 

02 

03 

0A 

05 

06 

07 

0S 

09 

010 


k  =  4 

0.3338537800974016E-01 
-0.1 7370662 12781 738E+00 
0.1198836762516400E+01 
-0.2591374963545573E+01 
0.3883697343902998E+01 
-0.2430220427598129E+01 
0.6270755476376530E+00 
-0.4769301964491491E-01 


ib  =  5 

0.1916264301337224E-02 
0.1567268574189051E+00 
-0.4334786391079948E+00 
0.7215818540402889E+00 
0. 1482966328206821 E+0 1 
-0.6595938685246877E+01 
0.1169258889232862E+02 
-0.9958531304701892E+01 
0.4103424553775595E+01 
-0.6712561210148048E+00 


Singularity  s(i)  =  (x  -  a) 
k  =  2  k  =  3 

01  0.8353164416920675E-01  0.3888639530245085E-01 

02  0.191 1747004918826E+00  0.8640672385266692E-01 

03  0.3657110693973134E+00  0.3880300543933000E-01 

04  -0.1404174140584028E+00  0.7069014304227469E+00 

0s  -0.45372631 16382184E+00 

00  0.8272875662102377E-01 
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Table  3.8:  (continued) 


A 

02 

03 

04 

05 

0e 

07 

0S 

09 

010 


Singularity  4(1)  =  (x  —  (continued) 

k=4  k=5 


0.2823063071751857E-01 

-0.1016030472473191E+00 

0.8957498450168363E+00 

-0.19446300il80514lE+01 

0.3079540564910779E+01 

-0.1835281366554477E+01 

0.3822511312525950E+00 

-0.4257746290791726E-02 


0. 160451 780304348 1 E  -0 1 
-0.2672763418512212E-01 
0.4952499866825667E+00 
-0.1916835699721412E+01 
0.6247645777407226E+0 1 
-0.12330931 15169777E+02 
0.1631924962150504E+02 
-0. 1237865585994330E+02 
0.4849695950284900E+0 1 
-0.7747361683625560E+00 


01 

02 

03 

04 

05 

00 


Singularity  ^(x)  =  (i  -  a)*/^ 
k=2  k=Z 


0.3685550676579945E-01 

0.3411680628277778E+00 

0.2026521530973060E+00 

-0.8067572269088319E-01 


0.1729631990952830E-01 
0.202444 1 1 14244492E+00 
-0.2155157764926719E+00 
0.9910326020047757E+00 
-0.6154646588009935E+00 
0.1202074019549121E+00 


=  4 

01  0.7512196681423a55E-01 

02  -0.4538891994788910E+00 

0z  0.2050275204711021E+01 

0^  -0.4084723895152133E+01 

A  0.5503263764388441E+01 

0s  -0.3512451256619964E+01 

07  0.1038747521068656E+01 

0s  -0.1163441057313696E+00 

09 
0\o 


k  =  5 

-0.6317627640571476E-01 
0.7370068868212822E+00 
-0.2828700163634487E+01 
0.6657003574957570E+01 
-0.8199869710972154E+01 
0.4165258132856930E+01 
0.3554455470298683E+01 
-0.5923097194990581  E+01 
0.2913478385611404E+01 
-0.5123591045429318E+00 


Chapter  4 

Vector  Space  Bases 


In  this  chapter  we  construct  a  class  of  bases,  analogous  to  the  multi-wavelet 
bases  of  Chapter  2,  that  transform  the  dense  n  x  n-matrices  resulting  from  the 
discretization  of  integral  equations  into  sparse  matrices  with  order  0(n  log  n)  non¬ 
zero  elements  (to  arbitrary  finite  precision).  In  these  bcises,  the  inverse  matrices 
are  also  sparse,  and  are  obtained  in  order  0(n  log^  n)  operations  by  the  classical 
Schulz  method. 

A  recent  paper  [3]  introduces  the  use  of  wavelets  for  the  application  of  an 
integral  operator  to  a  function  in  0(n  log  n)  operations,  where  n  is  the  number 
of  points  in  the  discretization  of  the  function.  In  the  present  chapter,  rather 
than  employ  a  wavelet  basis  for  we  construct  bases  in  which  we  represent 

the  operators  discretized  by  Nystrom’s  method.  If  5  =  denotes 

the  points  of  the  discretization,  we  define  a  cl2iss  of  wavelet-like  bases  for  the 
finite-dimensional  space  of  functions  defined  on  5. 


4.1  Wavelet  Bases 

4.1.1  Properties  of  the  Bases 

Given  a  set  of  n  distinct  points  5  =  {xi,X2,...,x„}  C  'R-  (the  discretization)  we 
construct  an  orthonormal  beisis  for  the  n-dimensional  space  of  functions  defined 
on  S.  For  simplicity,  we  assume  that  n  =  k-2‘,  where  k  and  /  are  natural  numbers, 
and  that  Xi  <  X2  <  •  •  •  <  Xn-  The  basis  heis  two  fundamental  properties: 

1.  All  but  k  basis  vectors  have  k  vanishing  moments,  and 

2.  The  basis  vectors  are  non-zero  on  different  scales. 

Fig.  4.1  illustrates  a  matrix  of  basis  vectors  for  n  =  128  and  A:  =  4.  Each  row 
represents  one  basis  vector,  with  the  dots  depicting  non-zero  elements.  The 
first  k  basis  vectors  are  non-zero  on  xi,...,X2it,  the  next  k  are  non-zero  on 
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X2k+i,  •  ■  -  ,X4kj  and  so  forth.  In  all,  one  half  of  the  basis  vectors  are  non-zero 
on  2k  points  from  5,  one  fourth  are  non-zero  on  4k  points,  one  eighth  are  non¬ 
zero  on  8k  points,  etc.  Each  of  these  n/2  -I-  n/4  -|-  n/8  +  ‘  ’  •  +  k  =  n  —  k  basis 
vectors  has  k  zero  moments,  i.e,,  if  6  =  (6i, . . . ,  6„)  is  one  of  these  vectors,  then 

n 

Xi^  =  Q,  j  =  0, 1, , . . ,  fc  -  1. 

i=l 

The  final  k  vectors  result  from  orthogonalization  of  the  moments  {xi^,X2^ , . . . ,  Xn^) 
for  j  =  0, 1, . . . ,  A:  —  1. 

These  properties  of  local  support  and  vanishing  moments  lead  to  efficient 
representation  of  functions  which  are  smooth  except  at  a  finite  set  of  singularities. 
The  projection  of  such  a  function  on  a  vector  of  this  basis  will  be  negligible 
unless  the  vector  is  non-zero  near  one  of  the  singularities.  As  a  simple  example, 
we  consider  the  function  f(x)  =  log(x)  on  the  interval  [0,  l]  with  the  uniform 
discretization  Xi  =  i In.  A  hand  calculation  shows  that  for  any  c  >  0,  /  may  be 
interpolated  on  the  interval  [c,  2c]  by  a  polynomial  of  degree  7  with  error  bounded 
by  4~®,  or  roughly  single  precision  accuracy.  If  we  choose  A:  =  8  in  constructing 
the  basis,  /  will  be  represented  to  this  accuracy  by  the  k  basis  vectors  non-zero 
on  ii, . . . , xjfc,  the  k  basis  vectors  non-zero  on  Xi, . . . , x^jk,  and  so  forth,  down  to 


Figure  4.1:  The  matrix  represents  a  wavelet  basis  for  a  discretization  with  128 
points,  for  k  =  4.  Each  row  denotes  one  basis  vector,  with  the  dots  depicting 
non-zero  elements.  All  but  the  final  k  rows  have  k  vanishing  moments. 
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the  k  basis  vectors  non-zero  on  in  addition  to  the  k  orthogonalized 

moment  vectors.  The  number  of  non-negligible  coefficients  in  the  expansion  of  / 
in  this  basis  grows  logarithmically  in  n,  the  number  of  points  of  the  discretization. 
Although  this  example  is  idealized,  its  behavior  is  representative  of  the  general 
behavior  of  an  analytic  function  near  a  singularity. 

4.1.2  Construction  of  the  Bases 

The  conditions  of  “local”  support  and  zero  moments  determine  the  basis  vectors 
uniquely  (up  to  sign)  if  we  require  somewhat  more  moments  to  vanish.  Namely, 
out  of  the  k  vectors  non-zero  on  Zi, . . .  ,Z2it,  we  require  that  one  have  k  vanishing 
moments,  a  second  have  A:  -I- 1,  a  third  have  k-\-2,  and  so  forth,  and  the  fcth  have 
2A:  —  1  vanishing  moments.  We  place  the  same  condition  on  the  k  basis  vectors 
non-zero  on  xjfc+i, . . . , Z4jt,  and  so  on,  for  each  block  of  k  beisis  vectors  among 
the  n  —  k  basis  vectors  with  zero  moments. 

We  construct  the  basis  by  construction  of  a  finite  sequence  of  bases  (shown  in 
Fig.  4.2),  each  obtained  by  a  number  of  orthogonalizations.  The  first  basis  results 
from  nf{2k)  Gram-Schmidt  orthogonalizations  of  2k  vectors  each.  In  particular, 
the  vectors  {xi^, . . . ,  x^V)  for  j  =  0, . . .  ,2A:  —  1  are  orthogonalized,  the  vectors 
(x2fe+i^, .  • . ,  Xik^)  for  j  =  0, . . . ,  2A;  —  1  are  orthogonalized,  and  so  forth,  up  to  the 
vectors  (x„_2fc+i-', . . . ,  Xn-')  for  j  =  0, . . . ,  2Ar  —  1  which  are  orthogonalized. 

Half  of  the  n  vectors  of  the  first  basis  have  at  least  k  zero  moments;  in 
forming  the  second  basis,  these  vectors  are  retained;  the  remaining  n/2  basis 
vectors  are  transformed  by  an  orthogonal  transformation  into  basis  vectors,  each 
of  which  is  non-zero  on  4k  of  the  points  zi,...,x„,  and  half  of  which  have  at 
least  k  vanishing  moments.  The  orthogonal  transformation  results  from  nl{4k) 
Gram-Schmidt  orthogonalizations  of  2k  vectors  each.  Similarly,  the  third  basis 
is  obtained  from  the  second  basis  by  an  orthogonal  transformation  that  itself 
results  from  n/{Sk)  Gram-Schmidt  orthogonalizations  of  2k  vectors  each.  Before 
we  can  specify  these  orthogonalizations,  we  require  some  additional  notation. 

Suppose  F  is  a  matrix  whose  columns  ui, . . . ,  V2k  are  linearly  independent.  We 
define  W=Orth(V)  to  be  the  matrix  which  results  from  the  column-by-column 
Gram-Schmidt  orthogonalization  of  V.  Namely,  denoting  the  columns  of  W  by 
we  have 

linear  spanfioj, . . . ,  lu,}  =  linear  span{i;j, . . . ,  u,}  •  •  _  i  ol 

For  a.  2k  X  2A:-matrix  V  we  let  and  denote  two  k  x  2A:-matrices, 
consisting  of  the  upper  k  rows  and  the  lower  k  rows  of  V, 
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Now  we  proceed  to  the  definition  of  the  basis  matrices.  Given  the  set  of  points 
S  =  {xi, . . . ,  x„}  C  'R  with  xi  <  •  •  •  <  Xn,  where  n  =  •  2^  we  define  the  2k  x  2k 

moments  matrices  Mi,,  for  i  =  1, . . . ,  nf{2k)  by  the  formula 


/  1  x,.+i  •••  x,.+i^*~^  \ 

1  x,,+2  •  •  • 

\  1  X*<+2it  •••  / 


(4.1) 


Figure  4.2:  Each  of  the  four  matrices  represents  one  basis,  as  in  Fig.  1.  The 
upper-left  matrix  is  formed  by  orthogonalizing  moment  vectors  on  blocks  of  2k 
points.  The  upper-right  matrix  is  obtained  from  the  upper-left  matrix  by  premul¬ 
tiplying  by  an  orthogonal  matrix  which  is  the  identity  on  the  upper  half.  Simi¬ 
larly,  the  lower  matrices  are  obtained  by  further  orthogonal  transformations.  The 
lower-right  matrix  represents  the  wavelet  basis  for  n  =  64,  k  =  4. 
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where  5,-  =  (i  —  l)2fc.  The  first  basis  matrix  Ui  is  given  by  the  formula 


Ui  = 


!  Ui.,'- 


TT  ^ 


Uxa 


U, 


l,ni 


V 


where  =  Orth(A/x,i)  and  ni  =  nf{2k).  The  second  basis  matrix  is  U^Ui, 
with  U2  defined  by  the  formula 


U, 


where  /„  is  the  m  x  m  identity  matrix  and  U2  is  given  by  the  formula 


t/;  = 


/  Ui.,‘ 


II  ^ 

^2,1 


II  ^ 

^2,2 


U2,n,^ 


where  n2  =  nl{4k),  U2/  =  Orth(M2,<),  and  M2,,  is  given  by 


M2..  = 


^I.2«-1^Mi,2,-1  \ 
^l,2i^Mi,2i  / 


In  general,  the  jth  basis  matrix,  for  j  =2,...,  log2(n/fc),  is  Uj  ■  •  -  Ui,  with  Uj 
defined  by  the  formula 


/n- 


n/2^-» 
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where  f/j  is  given  by  the  formula 

l 


11-  ^ 


Via’- 


^i.2 


V  ) 


where  rij  =  n/(2-’  A;),  f/j,,  is  given  by 

Uj/  =  Orth(Mj.<), 

and  A/j,,-  is  given  by 


M,, = ( ) . 


(4.2) 


(4.3) 


The  final  basis  matrix  U  =  Ur"  Ui,  where  /  =  log2(n/fc),  represents  the  wavelet 
basis  of  parameter  k  on  Xi,. ,  Xn. 


Remark  4.1  The  definitions  given  for  the  basis  matrices  are  mathematical  def¬ 
initions  only;  in  a  numerical  procedure,  considerable  roundoff  error  would  be 
introduced  by  the  orthogonalizations  defined  above.  In  the  actual  implementa¬ 
tion,  the  matrices  Mj,i  are  shifted  and  scaled,  resulting  in  a  numerically  stable 
procedure  that  is  equivalent  to  the  above  definitions  (in  exact  arithmetic).  Details 
of  this  procedure  are  provided  in  Section  4.3. 


It  is  apparent  that  the  application  of  the  matrix  U  to  an  arbitrary  vector 
of  length  n  may  be  accomplished  in  order  0(n)  operations  by  the  application 
of  Ui,. . .  ,Ui  in  turn.  Similarly,  U~^  =  may  be  applied  to  a  vector  in  order 
0(n)  operations.  Certain  dense  matrices,  in  particular  those  arising  from  integral 
operators,  are  sparse  in  the  basis  of  U  and  their  similarity  transformations  can 
be  computed  in  0(n  log  n)  operations.  These  techniques  are  developed  in  the 
following  sections. 

We  conclude  this  section  with  an  illustration  of  the  vectors  of  one  baisis  from 
this  class,  in  Fig.  4.3. 
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Figure  4.3:  Basis  vectors  on  four  scales  are  shown  for  the  basis  where  n  =  128, 
points  ii, . . .  ,z„  are  equispaced,  and  k  =  8.  The  first  column  of  vectors  consists 
of  rows  1-8  of  U ,  the  second  column  consists  of  rows  65-72,  etc.  Note  that  half 
of  the  vectors  are  odd  and  half  are  even  functions,  and  that  the  odd  ones  are 
generally  discontinuous  at  their  center. 


56 


CHAPTER  4.  VECTOR  SPACE  BASES 


4.2  Second-Kind  Integral  Equations 

4.2.1  Sparse  Representation  of  Integral  Operators 

As  in  Chapters  2  and  3,  we  concern  ourselves  here  with  kernels  K  =  K{x,  <)  which 
are  analytic  except  at  x  =  t,  where  they  may  possess  an  integrable  singularity. 
We  discretize  the  integral  operator  AC, 

by  Nystrom’s  method,  using  an  equispaced  quadrature.  Given  n  >  2,  we  define 
points  Xi, . . . ,  Xn  to  be  equispaced  on  the  interval  [a,  b], 

Xi  =  a  +  (i  -  1)(6- a)/(n  -  1),  (4.4) 

and  define  the  elements  T.-j  of  the  n  x  n-matrix  T  by  the  formula 

=  (4.5) 

Note  that  the  matrix  T  —  2  (n)  corresponds  to  a  primitive,  trapezoid- like  quadra¬ 
ture  discretization  of  the  integral  operator  AC.  The  matrix  T  possesses  the  same 
smoothness  properties  as  the  kernel  A'(x,<).  Transformation  of  T  by  the  bases 
developed  in  Section  4.1  produces  a  matrix  that  is  sparse,  to  high  precision.  The 
number  of  elements  is  effectively  bounded  by  order  O(nlogn). 

When  the  matrix  representing  the  quadrature  corrections  developed  in  Sec¬ 
tion  3.2  is  added  to  T,  producing  high-order  convergence  to  the  integral  operator, 
this  bound  remains  valid.  . 

The  matrix  T,  transformed  by  the  orthogonal  n  x  n-matrix  U,  can  be  decom¬ 
posed  into  the  sum  of  a  sparse  matrix  and  a  matrix  with  small  norm.  That  is, 
given  c  >  0,  there  exists  >  0,  independent  of  n,  such  that  the  transformed 
matrix  can  be  written  in  the  form 

UTU'^  =  V  +  E, 

where  the  number  of  elements  in  V  =  V{n)  is  bounded  by  Ct  n  log  n  and  E  =  E{n) 
is  small:  ||  Aj|  <  c  ||T||.  This  assertion  can  be  proven  very  similarly  to  Lemmas  2.3 
and  2.4;  in  fact,  substitution  of  the  finite  sums  which  determine  the  elements  of 
UTU^  for  the  integrals  in  those  lemmas  yields  the  assertion. 

4.2.2  Solution  via  Schulz  Method 

The  sparse  matrix  representing  the  integral  operator  also  has  a  sparse  inverse, 
as  was  the  case  for  the  multi-wavelet  bases  developed  in  Chapter  2.  We  invert 
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the  matrix  by  the  Schulz  method,  described  in  Section  2.4.3.  Each  iteration 
requires  two  matrix-matrix  multiplications,  where  the  multiplicands  are  sparse, 
and  yields  a  sparse  product.  The  method  is  quadratically  convergent;  the  number 
of  iterations,  given  by  bound  (2.23,  grows  only  as  the  logarithm  of  the  condition 
number. 


4.2.3  Oscillatory  Coefficients 

We  now  consider  a  somewhat  more  general  class  of  integral  equations,  in  which 
the  integral  operator  is  given  by  the  formula 

{DKf){x)  =  p{x)  I  K{x,  t)  f{t)  dt, 

where  the  kernel  K  is  assumed  to  be  smooth,  but  the  coefficient  function  p 
can  be  oscillatory.  In  particular,  we  only  restrict  p  to  be  positive.  In  terms  of 
generality,  these  problems  lie  between  the  problems  with  smooth  kernels  (and 
constant  coefficient)  and  those  with  arbitrary  oscillatory  kernels. 

Writing  the  corresponding  integral  equation  in  operator  form,  we  obtain  the 
equation 

{I-DK)f  =  g.  (4.6) 

Although  /}  is  a  diagonal  operator,  and  K  is  smooth,  it  is  clear  that  the  dis¬ 
cretization  of  the  operator  DK  will  not  be  a  sparse  matrix  in  wavelet  coordi¬ 
nates.  In  this  framework,  it  would  appear  that  the  constructions  of  Chapter  2 
and  this  chapter  are  inapplicable.  If  we  instead  consider  the  operator 
in  which  oscillations  in  the  rows  match  those  in  the  columns,  it  becomes  clear 
that  the  construction  of  Section  4.1  can  be  revised.  Rather  than  constructing 
basis  functions  orthogonal  to  low-order  polynomials  i-’,  we  can  construct  them 
to  be  orthogonal  to  p(i)‘/^i''.  The  sole  revision  in  our  definition  of  basis  matri¬ 
ces  I7i, ...  ,Ui  IS  to  replace  the  definition  (4.1)  of  the  moments  matrices  Mi,i  for 
i  =  1,. . .  ,n/(2A:),  by  the  new  definition 


/ 
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where  s,  =  (z  —  1)2^  and  pj  =  p(xjy^^. 

Now  the  integral  equation  (4.6)  can  be  transformed  to  the  equation 

(/  -  =  [D-^^^g). 

which  is  discretized  to  a  system  that  is  sparse  in  the  revised  wavelet  coordinates. 
The  inverse  matrix  is  also  sparse. 
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4.3  Numerical  Algorithms 

In  Section  4.1  we  defined  a  class  of  bases  for  functions  defined  on  {xi, . . . ,  x„}  and 
in  Section  4.2  we  showed  that,  to  finite  precision,  second-kind  integral  operators 
and  their  inverses  are  asymptotically  sparse  in  these  bases.  In  this  section  we 
present  procedures  for  computation  of  the  bases,  discretized  integral  operators 
in  these  bases,  and  the  inverses  of  these  operators.  In  Section  4.4  we  give  some 
numerical  examples  based  on  our  implementations  of  these  procedures. 

The  computation  of  the  wavelet  bases  is  discussed  next,  followed  by  a  dis¬ 
cussion  of  the  transformation  of  the  integral  operators  to  the  wavelet  bases.  We 
defer  discussion  of  the  computation  of  the  inverses,  sketched  above,  to  subsec¬ 
tion  4.3.3,  which  contains  detailed  descriptions  of  all  of  the  algorithms.  Finally, 
subsection  4.3.4  gives  the  complexity  analysis  for  the  algorithms. 

4.3.1  Computation  of  Wavelet  Bases 

It  was  mentioned  in  Section  4.1  that  the  mathematical  definition  of 
if  used  directly,  would  result  in  a  numerical  procedure  that  would  create  large 
roundoff  errors.  A  correct  procedure  is  obtained  by  shifting  and  scaling  the 
matrices  Mj,i  defined  there. 

For  a  pair  of  numbers  {ii,cr)  6  (7^\{0})  we  define  a2A:x2/:-matrix  5(/i,cr) 

whose  (t,i)th  element  is  the  binomial  term 

for  i  <  j,  and  S{n,a)ij  =  0  otherwise.  The  matrix  5(/i,cr)  is  upper-triangular 
and  non-singular,  and  its  inverse  is  given  by  the  formula 

=  S{-nl<7,  l/a).  (4.8) 

Furthermore,  the  product  formula 

5(/ii,cri)5(/X2,cr2)  =  S(/ii  -i-/i2<7i,cri<72)  (4.9) 

is  easily  verified. 

We  define  A/j  j  for  j  =  1, . . . ,  /  and  i  =  1, , . . ,  n/(2-'^)  by  the  formula 

M'.  =  A/,,S(^,,,<7,,),  (4.10) 

where  Hj,i  =  (xi+(,_i)fc2J  +Xifc2j)/2,  =  (x.wj  -  X\Mi-\)kv)l‘^  and  the  matrix 

is  defined  by  Eq.  (4.1)  and  Eq.  (4.3)  in  Section  4.1.  The  matrix  f/,,,  is  given 
by  the  formula 

U,/  =  Orth(A/',), 


(4.11) 
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which  is  equivalent  to  the  definition  given  by  Eq.  (4.2).  This  equivalence  imme¬ 
diately  follows  from  the  fact  that  S{fi,cr)  is  upper-triangular  and  non-singular. 
The  matrices  for  i  =  1, . . .  ,n/(2k)  are  actually  computed  by  the  formula 


Ki  = 


\  »l,i  / 

1  <ri,i  ) 

\  o\,i  ) 

(4.12) 


where  Sj  =  {i  —  1)2A:.  Likewise,  the  matrices  Mj,  for  j  =  2,...,/  and  i  = 
1, . . .  ,n/(2'’fc)  are  computed  by  the  formula 


"'■■“I  )' 


(4.13) 


where  5j  •  and  S?,  are  defined  by  the  formulae 

(4.14) 

(4.15) 

Application  of  the  inverse  and  product  rules  given  in  Eq.  (4.8)  and  Eq.  (4.9)  to 
Eq.  (4.14)  and  Eq.  (4.15)  yields  formulae  by  which  5],  and  Sj,  can  be  computed; 

-  fij.i,2i-i)l(Tj.i,2i.U  (rj,il<7j.i,2i-l)  (4.16) 

Sli  =  -  Hj.i,2i)l<rj.i,2i,  (4.17) 


The  matrices  Mj^  given  by  Eq.  (4.12)  and  Eq.  (4.13)  are  easily  seen  to  be 
mathematically  equivalent  to  those  defined  by  Eq.  (4.10);  nonetheless,  computa¬ 
tion  of  Mj,  using  Eqs.  (4.12)  and  (4.13)  avoids  the  large  roundoff  errors  which 
would  otherwise  result. 


4.3.2  Transformation  to  Wavelet  Bases 

We  assume  that  for  equispaced  points  zj, . . . ,  x„  (defined  in  Eq.  (4.4))  and  some 
k  the  orthogonal  matrices  f/i, . . . ,  C/)  defined  in  Section  4.1  have  been  computed 
(/  =  log2(n/A:)).  We  now  present  a  procedure  for  computation  of  UTU^,  where 
U  =  Ui  •  •  •  Ui  and  T  is  the  discretized  integral  operator  defined  in  Eq.  (4.5). 

Simple  Example 

We  begin  with  a  simplified  example  in  which  T  is  replaced  by  an  n  x  n-matrix 
V  of  rank  k  whose  elements  VJj  are  defined  by  the  equation 

rsl  3sl 


i,j  =  l,...,n. 
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Each  row  and  each  column  of  V  contain  elements  which  are  the  values  of  a 
polynomial  of  degree  k  —  1.  The  matrix  V  can  be  written  as  V  =  P^AP,  where 
the  elements  of  the  k  x  n-matrix  P  are  defined  by  Pij  =  Xj'~^  and  A  is  the 
k  X  fc-matrix  with  elements  A,j.  Recalling  that  the  last  k  rows  of  the  basis 
matrix  U  consist  of  an  orthogonalization  of  the  moment  vectors  , . . . ,  Xn^)  for 
j  =  0, . . . ,  A:  —  1,  we  can  rewrite  V  as  V  =  {P'Y A' P\  Here  the  k  x  n-matrix  P' 
consists  of  the  last  k  rows  of  U  and  A'  is  a  new  k  x  A:- matrix  with  elements  A\-. 

By  the  orthogonality  of  f/,  it  is  clea*  the  nx n-matrix  UVU^  =  U{P')^ A' P'U^ 
consists  entirely  of  zero  elements  except  the  k  x  A:-submatrix  in  the  lower-right 
corner,  which  is  the  matrix  A'.  Given  a  function  to  compute  elements  of  the 
n  X  n-matrix  V,  the  matrix  A'  can  be  computed  in  time  independent  of  n  by 
using  a.  k  X  k  extract  of  values  from  V.  We  form  the  matrix  V  with  elements  V^' 
defined  by  the  formula 

y'ij  =  yinikjnik-,  hj  =  1, . . . ,  A:,  (4-18) 

Then  V  =  {P")^ A' P",  where  P"  is  the  k  x  k  extract  of  P'  with  elements  given 
—  ^!,jn/k-  Thus  we  obtain 


A'  =  ((Py)->V"(P")'' 


(4.19) 


from  P"  and  V'  readily  in  0{k^)  operations,  and  we  have  obtained  UVU^ . 


General  Case 

The  integral  operator  matrix  T  is,  of  course,  not  of  low  rank,  but  it  can  be  divided 
into  submatrices,  each  approximately  of  rank  k  (see  Fig.  4.4).  The  submatrices 
near  the  main  diagonal  are  of  size  k  x  k,  those  next  removed  are  2k  x  2k,  and 
so  forth  up  to  the  largest  submatrices,  of  size  n/4  x  n/4.  The  total  number  of 
submatrices  is  proportional  to  n/k.  Given  an  error  tolerance  €  >  0,  A:  may  be 
chosen  (independently  of  n)  so  that  each  submatrix  of  T,  say  T\  may  be  written 
ais  a  sum,  T'  =  V'  +  E’,  where  the  elements  of  V  are  given  by  a  polynomial  of 
degree  A:  —  1  and  |1P'1|  <  cl|7’’||. 

The  simplified  example,  in  which  the  matrix  to  be  transformed  is  of  rank 
k,  is  now  applicable.  Each  submatrix  of  T  is  treated  as  a  matrix  of  rank  k 
and  is  transformed  to  wavelet  coordinates  (for  its  own  scale)  in  order  0(A:^) 
operations.  To  make  this  precise,  we  write  T  =  To  +  •  •  •  +  Ti-i  where  P,  consists 
of  the  submatrices  of  size  2‘A:  x  2*k.  For  each  i,  the  submatrices  of  T,  may  be 
interpolated  by  rank  k  submatrices,  as  indicated  by  the  extract  of  Eq.  (4.18),  to 
obtain  matrices  K.  Thus  Ti  =  Vi  +  Ei,  where  ||P,||  is  small.  In  the  simplified 
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Figure  4.4:  The  matrix  represents  a  discretized  integral  operator  with  a  kernel 
that  is  singular  along  the  diagonal.  The  matrix  is  divided  into  submatrices  of 
rank  k  (to  high  precision)  and  transformed  to  a  sparse  matrix  with  O(nlogn) 
elements.  Here  n/k  =  32. 

example  above,  we  have  shown  that  the  transformed  matrices 


VFo 

=  K) 

W, 

= 

W2 

(4.20) 

W,.2 

=  Ul.2---U,Vi.2Ux^  ■■■Ul.2'^ 

can  be  computed  by  many  applications  of  Eq.  (4.19),  all  in  order  Olnk^)  oper¬ 
ations.  This  estimate  follows  from  the  fact  that  there  are  0{n/k)  submatrices, 
each  of  which  is  transformed  in  0{k^)  operations.  Now  we  define  n  x  n-matrices 
Ro, . .  .,Ri  recursively: 


\  Wo  i  =  0 

[UiR,.^U,'^  +  W,  z>l 


(4.21) 


(here  VF/_i  =  Wi  =  0).  Then  Rt  contains  the  final  result,  Ri  =  U{T  —  E)U^ , 
where  E  =  Eq  +  •  ■  ■  ■¥  Ei-2- 
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The  matrix- matrix  products  in  the  definition  of  /2o, can  be  computed 
directly,  since  the  factors  and  the  products  contain  no  more  than  O(nlogn)  ele¬ 
ments.  A  simple  implementation  with  standard  sparse  matrix  structures  results 
in  a  total  operation  count  of  order  O(nlog^n),  but  an  implementation  using 
somewhat  more  elaborate  data  structures,  in  which  repetitive  handling  of  data 
is  avoided,  requires  only  order  O(nlogn)  operations. 

Computation  using  the  result  Ri  is  made  more  efficient  by  removing  the  ele¬ 
ments  of  Ri  which  can  be  neglected,  within  the  precision  with  which  /?/  approx¬ 
imates  UTU^.  For  a  given  precision  e,  we  discard  a  matrix  E'  by  eliminating 
elements  from  Ri  below  a  threshold  r.  The  threshold  depends  on  the  choice  of 
norm;  in  our  implementation,  we  use  the  row-sum  norm 

II  All  =  max^lA.jl, 

for  an  n  X  n-matrix  A.  The  element  threshold 

r  =  -i||T|l  (4.22) 

n 

clearly  results  in  a  discarded  matrix  E'  with  ||£'||  <  e||T||. 

4.3.3  Detailed  Descriptions  of  Algorithms 

Procedure  to  compute  Ui,..  .,Ui 

Comment  [Input  to  this  procedure  consists  of  the  number  of  points  n, 
the  number  of  zero  moments  k,  and  the  points  Xj, . . . ,  x„.  Output  is  the 
matrices  t/,,,  for  j  =  1, . . . ,  /  and  1  =  1,...  ,n/(2-'/:),  which  make  up  the 
matrices  Ui, . . .  ,Ui  (note  I  =  log2(n/^:)).j 

Step  1. 

Compute  the  shifted  and  scaled  moments  matrices  AI[ ,  for 
i  =  1, . . .  ,n/(2A:)  according  to  Eq.  (4.12). 

Step  2. 

Compute  f/i,,  from  M[-  by  Eq.  (4.11)  using  Gram-Schmidt 
orthogonalization  for  i  =  .  ,nf{2k). 

Step  3. 

Comment  [Compute  M'-  and  Uj^,  for  j  —  2, . . .  J  and 
z  =  l,...,n/(2U-).] 
do  j  =  2, . . . ,  / 

do  z  =  1, . . .  ,n/{2^k) 
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Compute  M'j_i  2i-i  and  jj- 

Compute  5],-  by  Eq.  (4.16)  and  by  Eq.  (4.17); 

multiply  to  obtain  by  Eq.  (4.13). 
Orthogonalize  Mj  ,-  to  obtain  by  Eq.  (4.11). 
enddo 
enddo 


Procedure  to  compute  UTU^ 

Comment  [Input  to  this  procedure  consists  of  n,  k,  the  matrices  Uj^i 
computed  above,  a  function  to  compute  elements  of  T,  and  the  chosen 
precision  e.  Output  is  a  matrix  Ri  such  that  ||i?/  —  UTU'^\\  <  e||r||.] 

Step  4. 

Compute  the  k  x  k  extracts,  indicated  by  Eq.  (4.18),  of  the 
submatrices  of  T  shown  in  Fig.  4.4. 

Step  5. 

Extract  the  matrices  P"  (Eq.  (4.19))  from  f/j,  U2U1,  . . . ,  Ut  •  ■  'Ui 
and  compute  Wq,  . . . ,  according  to  Eqs.  (4.20). 

Step  6. 

Compute  Rq,  . . .  ,Ri  by  Eq.  (4.21),  discarding  elements  below 
a  threshold  r  determined  by  the  precision  c  (Eq.  (4.22)). 


Procedure  to  compute  UT~^U^ 

Comment  [Input  to  this  procedure  consists  of  n,  the  matrix  Ri  which 
approximates  UTU^  ,  and  the  precision  e.  Output  is  a  matrix  Xm 
that  approximates  UT~^U^.] 


Step  7. 

Compute  the  matrix  Xq  =  Ri^ Ri/\\Rt^ Ri\\  by  direct  matrix 
multiplication,  discarding  elements  below  a  threshold  r 
determined  by  the  precision  e  (Eq.  (4.22)). 

Step  8. 

Comment  [Obtain  the  inverse  by  Schulz  iteration.] 
do  m  =  0, 1, . . .  while  ||/  -  >  c 

Compute  Xm+i  =  2Xm  —  X^RiX^,  discarding  elements 
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below  threshold, 
enddo 


4.3.4  Complexity  Analysis 

In  the  following  table,  we  provide  the  operation  count  for  each  step  of  the  com¬ 
putation  of  UT~^U^. 

Step  Complexity  Explanation _ 

1.  0{nk)  There  are  n/{2k)  2k  x  2A:-matrices;  each  element 

of  the  matrices  is  computed  in  constant  time. 

2.  0{nk^)  For  each  of  the  nf{2k)  matrices,  perform  a  Gram- 

Schmidt  orthogonalization  requiring  order  O(k^) 
operations. 

3.  0(nA:^)  For  each  of  n/ (4k) +  n/ (8k) -I - hl=n/(2k)  —  l 

matrices,  compute  four  products  of  a  k  x  2k- 
matrix  with  a  2k  x  2l:-matrix,  construct  two 
2k  X  2fc-matrices,  and  orthogonalize  one  2k  x  2k- 
matrix. 

4.  0(nk)  There  are  6(1+3+74 - |-(n/(2A:)  — 1))+ 3(n/fc)  — 

2,  or  order  0(n/k),  submatrices  of  T  and  for  each 
matrix  we  compute  P  elements. 

5.  0(nk^)  There  are  nf  (2k)  4-71  f  (4k)  4 - t-l=n/fc  —  l  ma¬ 

trices  /*",  each  the  product  of  two  k  x  Ar-matrices. 
These  are  each  inverted  and  multiplied  with  the 
0(nfk)  matrices  of  the  previous  step. 

6.  0(n  log  n)  The  diagonally-banded  matrix  Wq,  which  contains 

0(n)  elements,  grows  to  O(nlogn)  elements  by 
the  computation  of  UWqU^ ,  eis  can  be  seen  by 
simply  examining  pictures  of  Wq  and  U.  The  non¬ 
zero  elements  of  the  transformed  l-Fj , . . . ,  are 
a  subset  of  those  of  Wq. 

7.  O(nlog^n)  Multiplication  of  two  matrices,  each  with  order 

O(nlogn)  elements,  to  obtain  a  product  with  or¬ 
der  O(nlogn)  elements. 
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Step  Complexity  Explanation _ 

si  0(n  log'*  n)  Two  multiplications  like  that  of  Step  7  are  made 
per  iteration;  the  number  of  iterations  is  indepen¬ 
dent  of  n  and  given  by  bound  (2.23). 

Total  0(r  log^  n) 

4.4  Numerical  Examples 

In  this  section  we  present  operators  from  several  integral  equations,  the  discretiza¬ 
tion  and  transformation  of  the  operators  to  our  wavelet  bases,  and  the  inversion 
of  the  operators  via  Schulz  method. 

4.4.1  Uncorrected  Quadratures 

We  first  examine  simple  quadratures  with  equal  weights,  except  weight  zero  at 
the  singularity,  as  represented  by  matrix  T  =  T{n)  defined  by  Eq.  (4.5).  We 
transform  the  matrix  /  —  T  to  wavelet  coordinates  as  described  in  Section  4.3.2, 
then  compute  (/  —  T)~^. 

These  discretizations  are  not  particularly  useful  for  the  solution  of  the  inte¬ 
gral  equations,  due  to  their  slow  convergence  to  the  integral  operators.  They 
nonetheless  make  good  illustrative  examples,  for  they  retain  the  smoothness  of 
the  operator  kernels  and  produce  correspondingly  sparse  matrices.  In  the  next 
subsection,  we  examine  the  results  of  using  high-order  quadratures. 

For  various  sizes  n  of  discretization,  we  tabulate  the  average  number  of  ele¬ 
ments  per  row  in  the  transformed  matrix  U{I  —  T)U^  and  the  computation  time 
to  obtain  the  matrix.  In  addition,  we  display  the  average  number  of  elements 
per  row  of  its  inverse,  and  the  time  to  compute  the  inverse.  Finally,  we  show  the 
error  introduced  by  these  computations.  The  error  is  determined  by  the  appli¬ 
cation  of  the  forward  and  inverse  transformations  to  a  random  vector:  Choose  a 
vector  V  of  length  n  with  uniformly  distributed  pseudo-random  elements;  com¬ 
pute  {I  —  T)v  directly,  by  a  standard  procedure  requiring  order  O(n^)  operations; 
transform  to  wavelet  coordinates,  obtaining  U(I  —  T)v,  apply  the  computed  value 
of  U{I  —  T)~^U^  to  the  vector  U{I  —  T)v\  transform  to  original  coordinates  by 
application  of  f/^;  compare  the  result  to  v.  The  measure  of  error  is  the  relative 

error  c/;j,  defined  by  Eq.  (3.27). 

The  programs  to  transform  and  invert,  2is  well  as  those  to  determine  the  error, 
were  implemented  in  FORTRAN.  All  computations  were  performed  in  double¬ 
precision  arithmetic  on  a  Sun  Sparcstation  1. 

The  first  set  of  examples  is  for  the  kernel  K[x,  t)  =  log  [x  —  t\,  for  a  wavelet 
basis  of  order  =  4  and  various  choices  of  precision  e.  The  matrix  sparsities, 
execution  times,  and  errors  appear  in  Table  4.1.  Although  the  sparse  matrices 
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Table  4.1:  The  integral  operator  JC  defined  by  the  formula  [Kf){x)  =  f{x)  — 
/o  log  |x  —  t\f{t)dt  is  discretized,  transformed  to  the  wavelet  coordinates  with 
k  =  4,  and  inverted.  For  various  precisions  e  and  various  sizes  of  discretization, 
we  tabulate  the  average  number  of  elements/row  Ni  of  the  matrix  in  wavelet 
coordinates  and  the  time  in  seconds  ti  to  compute  it,  corresponding  statistics  N2 
and  <2  for  the  inverse,  and  the  error  (see  text). 


Transform.  Inversion 


£2 


€ 

n 

Nx 

ti 

N2 

h 

Error 

lo-’-* 

64 

7.2 

2 

8.3 

2 

0.503E-02 

128 

5.9 

3 

6.5 

4 

0.257E-02 

256 

3.8 

7 

4.4 

4 

0.250E-02 

512 

2.8 

13 

3.1 

6 

0.236E-02 

1024 

1.9 

26 

2.1 

6 

0.227E-02 

2048 

1.4 

49 

1.4 

6 

0.221E-02 

4096 

1.2 

97 

1.2 

8 

0.221E-02 

8192 

1.1 

195 

1.1 

12 

0.217E-02 

10-3 

64 

17.6 

2 

19.5 

14 

0.350E-03 

128 

18.1 

5 

20.0 

36 

0.270E-03 

256 

18.0 

11 

20.0 

83 

0.331E-03 

512 

14.5 

21 

15.7 

123 

0.257E-03 

1024 

13.3 

41 

15.5 

262 

0.340E-03 

2048 

8.5 

73 

9.8 

287 

0.233E-03 

4096 

5.8 

131 

6.5 

304 

0.222E-03 

8192 

3.7 

242 

4.4 

312 

0.221E-03 

10-^ 

64 

28.4 

3 

30.3 

36 

0.104E-03 

128 

32.1 

6 

34.3 

111 

0.140E-03 

256 

34.5 

15 

37.5 

302 

0.161E-03 

512 

33.1 

31 

35.8 

618 

0.177E-03 

1024 

30.2 

63 

33.6 

1280 

0.189E-03 

2048 

25.0 

121 

27.6 

2040 

0.192E-03 

banded,  we  loosely  refer 

to  the 

average 

number 

of  matrix  elements  per 

row  as  the  matrix  bandwidth.  We  make  the  following  observations; 

1.  The  bandwidths  Ni,  N2  of  the  operator  and  its  inverse  decrease  with  increas¬ 
ing  matrix  size.  In  other  words,  in  the  range  of  matrix  sizes  tabulated,  the 
number  of  matrix  elements  grows  sublinearly  in  the  matrix  dimension  n. 

2.  The  operator  matrix  in  wavelet  coordinates  is  computed  in  time  that  grows 
nearly  linearly  in  n. 
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Figure  4.5:  The  matrices  constructed  in  the  transformation  of  I  —  T,  matrices 
Ro,...,R3  defined  in  Eq.  (4.21),  are  shown  for  kernel  K{x,t)  =  log  |x  —  f|, 
e  =  10~^,  and  n  =  64.  The  matrix  R4  looks  like  Rj  and  is  not  shown. 

3.  The  inverse  matrix  is  computed  in  time  which  grows  sublinearly  in  n.  This 
is  due  to  the  fact  that  the  cost  of  multiplying  the  sparse  matrices  is  roughly 
order  O(niV^),  for  size  n  and  bandwidth  N.  One  result  is  that  the  cost 
sometimes  drops  cis  n  incre«ises. 

4.  The  accuracy  is  within  the  precision  specified.  In  fact,  due  to  the  conser¬ 
vative  element  thresholding  (Eq.  4.22),  the  actual  error  is  considerably  less 
than  c. 

5.  The  cost  increases  with  increasing  precision  e,  due  to  the  increasing  band- 
widths  generated.  The  the  bandwidths  incretise  approximately  as  log(l/e). 

6.  For  k  =  4,  our  fast  transformation  algorithm  does  not  maintain  the  specified 
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Figure  4.6:  Transformed  matrix  U{I  —  T)U^  (^op)  and  its  inverse  (bottom)  are 
shown  for  kernel  K{x,t)  =  log  \x  —  e  =  10"^,  and  n  =  128. 


precision  of  c  =  10  This  anticipated  result  follows  from  the  error  estimate 
for  polynomial  interpolation  of  logarithm  on  intervals  separated  from  the 
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origin.  An  unanticipated  attendant  result  is  that  the  bandwidth  increases 
as  the  quality  of  approximation  deteriorates  (compare  to  fc  =  8,  below).  As 
a  result,  we  did  not  complete  examples  for  n  =  4096,  8192. 

7.  The  inversion  of  the  8192  x  8192  matrix  preserving  3-digit  accuracy  is  done 
in  5  minutes  on  the  Sparcstation.  This  compares  to  95  days  (estimated) 
for  inverting  the  dense  matrix  by  Gauss-Jordan  and  to  24  minutes  for  one 
dense  matrix- vector  multiplication  of  that  size. 

The  condition  number  of  the  problem,  as  approximated  by  the  product  of  the 
row-sum  norms  of  U{I  —  T)U^  and  its  computed  inverse,  is  3  (independent  of 
size).  Five  iterations  were  required  by  the  Schulz  method  to  achieve  convergence. 

In  Fig.  4.5  we  show  stages  in  the  transformation  of  the  matrix  I  —  T.  In  par¬ 
ticular,  for  c  =  10~^  and  n  =  64,  the  matrices  Rq^  . . . ,  defined  in  Eq.  (4.21) 
are  shown.  In  addition,  for  n  =  128  the  transformed  matrix  U{I  —  T)U^  and  its 
inverse  are  shown  in  Fig.  4.6. 

In  the  next  set  of  examples,  for  which  results  are  displayed  in  Table  4.2,  we 
used  the  wavelet  basis  of  order  k  =  S.  We  observe: 

1.  The  bandwidths  of  the  operator  matrix  and  its  inverse  are  less  for  fc  =  8 
than  for  A:  =  4.  The  inversion  times  are  correspondingly  smaller. 

2.  The  time  required  to  compute  the  operator  matrix  is  almost  four  times 
as  large  as  that  for  k  —  4.  This  is  due  to  the  cost  of  transforming  the 
near-diagonal  band,  which  is  twice  as  wide  for  A:  =  8  as  for  A:  =  4. 

3.  The  obtained  accuracy  exceeds  the  specified  precision  consistently. 

4.  As  for  k  =  4,  the  scaling  with  size  n  is  linear  for  the  transformation  step 
and  sublinear  for  the  inversion  step. 

In  the  final  set  of  examples  in  which  uncorrected  quadratures  were  used,  we 
perform  computations  for  A:  =  4  and  e  =  10“^,  with  various  operator  kernels. 
Table  4.3  presents  the  results.  The  first  three  kernels  contain  singularities  of  the 
types  s(i)  =  log(x)  and  s(x)  =  x°  for  a  =  ±5,  and  are  nonsymmetric  and  non- 
convolutional.  It  is  readily  seen  that  the  bandwidth  is  strongly  dependent  on  the 
type  of  singularity,  with  the  singularity  x“’/^  producing  the  greatest  bandwidth. 
We  mention  also  that  this  particular  integral  equation  is  poorly-conditioned; 
the  condition  numbers  of  the  discretizations  for  n  =  64,128,256,512,1024  are 
9,17,34,98,469,  respectively. 

The  fourth  kernel  provides  an  example  with  an  oscillatory  coefficient  p(x)  = 
(1  +  |sin(100x)).  The  bases  developed  in  Section  4.2.3,  which  depend  on  p,  are 
used  to  transform  the  discretized  integral  operator  to  sparse  form.  We  see  in 
Table  4.3  that  the  inverse  is  also  very  sparse. 
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Table  4.2:  The  integral  operator  K  defined  by  the  formula  {K f){x)  =  f{x)  — 
/o  log  |x  —  t\f{t)dt  is  discretized,  transformed  to  the  wavelet  coordinates  with 
k  =  8,  and  inverted.  (See  Table  4^4  and  text.) 


Transform.  Inversion 


£2 


=T- 


n 

h 

N2 

<2 

Error 

64 

5.8 

4 

6.2 

1 

0.191E-02 

128 

5.0 

10 

5.-5 

2 

0.368E-02 

256 

3.3 

22 

3.e 

3 

0.184E-02 

512 

2.7 

46 

2.9 

4 

0.113E-02 

1024 

1.8 

92 

1.8 

4 

0.177E-02 

2048 

1.4 

182 

1.4 

5 

0.170E-02 

4096 

1.2 

363 

1.2 

8 

0.928E-03 

8192 

1.1 

729 

1.1 

11 

0.166E-02 

64 

13.4 

5 

14.5 

8 

0.373E-03 

128 

14.2 

13 

15.5 

21 

0.332E-03 

256 

13.5 

28 

14.5 

46 

0.259E-03 

512 

12.7 

57 

13.6 

90 

0.225E-03 

1024 

10.2 

114 

11.1 

134 

0.198E-03 

2048 

7.7 

221 

8.3 

176 

0.179E-03 

4096 

4.9 

429 

5.2 

185 

0.174E-03 

8192 

3.5 

818 

3.7 

208 

0.173E-03 

64 

21.8 

6 

23.0 

23 

0.280E-04 

128 

26.3 

15 

28.0 

81 

0.253E-04 

256 

28.7 

35 

31.0 

235 

0.246E-04 

512 

28.4 

75 

30.9 

538 

0.184E-04 

1024 

25.5 

149 

27.2 

969 

0.925E-05 

2048 

22.0 

297 

23.8 

1739 

0.899E-05 

4096 

17.7 

561 

19.1 

2610 

0.798E-05 

10 


10 


-3 


10- 


4.4.2  Solution  of  Integral  Equations 

In  the  preceding  subsection,  we  examined  the  characteristics  of  various  integral 
operators  and  their  inverses  in  wavelet  coordinates.  We  used  completely  straight¬ 
forward  discretizations;  the  quadratures  represented  sums  of  the  integrands  at 
equispaced  points  (excluding  singular  points).  Such  simple  quadratures  converge 
too  slowly  to  the  integral  operators  to  be  of  much  use  in  solving  integral  equa¬ 
tions,  and  we  now  turn  to  the  high-order  quadratures  developed  in  Chapter  3. 

We  first  present  examples  which  correspond  to  the  various  kernels  already 
tested  and  shown  in  Table  4.3.  In  Table  4.4  we  tabulate  the  results,  and  differ- 
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Table  4.3:  The  integral  operator  K.  defined  by  the  formula  (/C/)(x)  =  f{x)  — 
/o  K{x,t)  f{t)dt,  for  nonsymmetric,  nonconvolutional  kernels  K{x,t)  shown  be¬ 
low,  is  discretized,  transformed  to  the  wavelet  coordinates  with  k  =  4  and 


cos(if2)|x  _  64  27.2 

128  31.6 

256  35.6 

512  37.3 

1024  34.5 

cos(ii^)|x  -  64  6.8 

128  4.4 

256  2.9 

512  2.1 

1024  1.5 

2048  1.4 

4096  1.1 

8192  1.1 

(1  +  ^sin(100i))x  64  30.5 

Iog|x-f(  128  31.8 

256  21.2 

512  18.6 

1024  15.8 

2048  10.6 

4096  6.4 

8192  4.0 


3 

28.9 

32 

0.256E-03 

7 

34.1 

122 

0.357E-03 

16 

40.6 

454 

0.434E-03 

35 

46.3 

1509 

0.643E-03 

72 

45.4 

4166 

0.821E-03 

2 

7.3 

2 

0  303E-03 

4 

4.7 

2 

0.204E-03 

8 

3.0 

3 

0.209E-03 

15 

2.3 

3 

0.165E-03 

30 

1.5 

3 

0.208E-03 

60 

1.4 

6 

0.909E-03 

119 

1.2 

7 

0.614E-03 

242 

1.1 

12 

0.666E-03 

3 

33.8 

44 

0.344E-03 

6 

35.1 

103 

0.363E-03 

12 

24.1 

119 

0.348E-03 

23 

20.7 

225 

0.372E-03 

45 

18.4 

404 

0.392E-03 

82 

12.2 

466 

0.355E-03 

145 

7.4 

497 

0.336E-03 

265 

4.6 

510 

0.331E-03 

ences  from  Table  4.3  reflect  the  effect  of  the  quadratures. 

For  the  remaining  examples  we  choose  integral  equations  that  can  be  solved 
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Table  4.4:  The  integral  operator  K  defined  by  the  formula  (/C/)(x)  =  f{x)  — 
/q  K(x,t)  /(t)  dt,  for  nonsymmetric,  nonconvolutional  kernels  K{x,t).  shovm  be¬ 
low,  is  discretized  with  the  corrected  trapezoidal  rules,  transformed  to  the  wavelet 
coordinates  with  k  =  4  and  e  =  10~^,  and  inverted.  (Compare  to  Table  4-3.) 

Transform.  Inversion 

K{x,t)  n  Ni  ti  N2  t2  Error 

cos(it'')logli-  ^1  64  2O  4  sTs  38  0.164E-03 

128  31.5  9  34.3  103  0.162E-03 

256  30.8  21  33.9  221  0.172E-03 

512  27.0  41  29.7  370  0.177E-03 

1024  21.0  80  23.7  454  0.357E-03 

2048  14.8  143  17.2  566  0.317E-03 

4096  9.5  250  10.4  555  0.282E-03 

8192  5.8  448  6.9  665  0.271E-03 

cos(i<2)|i  -  <|-»/2  64  32.4  4  39.8  87  0.133E-02 

128  38.3  10  45.7  251  0.412E-03 

256  42.7  23  49.3  638  0.464E-03 

512  45.1  51  51.3  1494  0.562E-03 

1024  46.2  no  52.1  3309  0.635E-03 

cos(xf2)|x  -  64  10.4  3  18.4  9  0.867E-03 

128  7.6  6  13.8  13  0.526E-03 

256  5.1  13  9.3  16  0.358E-03 

512  3.3  25  5.2  15  0.292E-03 

1024  2.3  48  3.1  15  0.201E-03 

2048  1.9  96  2.3  20  0.393E-03 

4096  1.5  188  1.7  25  0.405E-03 

8192  1.3  374  1.4  36  0.404E-03 


analytically,  so  that  the  accuracy  of  the  method  can  be  checked.  We  consider  a 
class  of  integral  equations  with  logarithmic  kernel, 

\og  \x  -  t\f{t)dt  =  gm{x),  xG[0,  1],  (4.23) 

Jo 

where  the  right  hand  side  gm  is  chosen  so  that  the  solution  /  is  given  by  the 
formula  /(i)  =  sin(mx).  The  integration  can  be  performed  explicitly,  yielding 
the  formula 
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Table  4.5:  The  integral  equations  /(x)— log  |x  — f{t)  dt  =  for  which  an 

explicit  solution  is  known,  are  solved  by  the  methods  of  this  chapter  (compare  to 
Table  4-1  see  text).  For  e  =  10“^,  10~^,  10“'*  we  set  k  =  4,4,8,  respectively. 


Transform.  Inversion 


10-3 


128 

10.7 

7 

13.2 

14 

0.212E-02 

236 

8.6 

13 

10.6 

20 

0.140E-02 

512 

6.3 

26 

7.6 

26 

0.112E-02 

1024 

3.6 

48 

4.5 

28 

0.821E-03 

2048 

1.9 

90 

2.3 

21 

0.932E-03 

4096 

1.3 

174 

1.5 

15 

0.674E-03 

8192 

1.1 

344 

1.1 

13 

0.499E-03 

64 

27.7 

4 

31.3 

36 

0.235E-03 

128 

31.0 

9 

34.2 

99 

0.169E-03 

256 

30.6 

20 

33.6 

215 

0.161E-03 

512 

27.5 

41 

30.2 

377 

0.130E-03 

1024 

21.7 

79 

24.4 

470 

0.597E-03 

2048 

15.5 

143 

18.1 

604 

0.479E-03 

4096 

9.7 

248 

10.6 

579 

0.415E-03 

8192 

6.0 

444 

7.3 

690 

0.354E-03 

64 

37.2 

8 

45.9 

78 

0.127E-03 

128 

47.1 

23 

56.5 

278 

0.473E-04 

256 

52.9 

54 

60.9 

745 

0.311E-04 

512 

55.0 

118 

61.4 

1701 

O.lOOE-04 

1024 

52.3 

248 

57.2 

3287 

0.734E-05 

/  log  |x  —  <1  sin(mf)  =  log(x)  —  cos(m)  log(l  —  i) 

Jo 

—  cos(mx)[Ci(mx)  —  Ci(m(l  —  x))] 

—  sin(mx)[Si(mx)  +  Si(m(l  —  x))], 

where  Ci  and  Si  are  the  cosine  integral  and  sine  integral  (see,  e.g.,  [1],  p.  231). 
Equation  (4.23)  clearly  requires  quadratures  with  increasing  resolution  as  m  in¬ 
creases;  for  our  examples  we  let  n  —  m,  which  corresponds  to  27r  points  per 
oscillation  of  the  right  hand  side  gm- 

Initially  we  choose  coefficient  p(x)  =  1.  The  results  are  given  in  Table  4.6. 
Here  the  error  shown  is  the  error  of  the  computed  solution  relative  to  the  true 
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solution  of  the  integral  equation.  Many  of  the  observations  of  the  preceding 
examples  can  be  repeated  here;  additionally,  we  make  the  following  comments: 

1.  The  bandwidths  are  greater  than  for  the  uncorrected  quadratures,  but  this 
effect  generally  decreases  with  increasing  size. 

2.  The  integral  equations  are  solved  to  within  the  specified  precision  in  every 
case  but  one.  The  exception,  for  c  =  lO""*  and  n  =  64,  is  likely  due  to  the 
small  number  of  quadrature  points  and  high  specified  precision. 

3.  An  integral  equation  requiring  an  8192-point  discretization  is  solved  to  3- 
digit  accuracy  in  less  than  20  minutes  on  the  Sparcstation. 

For  our  second  set  of  integral  equations,  we  let  the  coefficient  p  be  the  oscil¬ 
latory  function  given  by  the  formula  p{x)  =  1  -f  5sin(100x).  We  carry  out  the 
transformation  described  in  Section  4.2.3  to  solve  the  integral  equation  4.23.  The 
results  are  shown  in  Table  4.6  and  as  with  Table  4.5  the  error  refers  to  the  error  of 
the  computed  solution  relative  to  the  true  solution  of  the  integral  equation.  For 
the  oscillatory  coefficient  we  see  performance  similar  to  the  constant-coefficient 
problem,  but  the  cost  is  higher. 
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Table  4.6:  The  integral  equations  f(x)  —  p{x)  /J  log  |x  —  f{t)  dt  =  gm{^),  for 
which  an  explicit  solution  is  known,  are  solved  by  the  methods  of  this  chapter 
(compare  to  Table  4-1  ond  see  text).  For  e  =  10“^,  10~^,  10“'*  we  set  k  =  4,4,3, 
respectively. 


10-3 


10-* 


Transform.  Inversion 


£2 


ji,  m 

Ni 

Ni 

<2 

Error 

64 

19.7 

4 

23.9 

18 

0.360E-02 

128 

17.7 

8 

21.0 

36 

0.182E-02 

256 

12.6 

15 

14.6 

47 

0.174E-02 

512 

8.4 

29 

9.8 

57 

0.112E-02 

1024 

4.7 

55 

5.7 

56 

0.104E-02 

2048 

2.4 

103 

2.7 

45 

0.902E-03 

4096 

1.6 

198 

1.7 

38 

0.720E-03 

8192 

1.3 

392 

1.3 

35 

0.543E-03 

64 

36.2 

4 

41.3 

63 

0.228E-02 

128 

40.8 

10 

47.0 

186 

0.209E-03 

256 

40.5 

23 

47.3 

427 

0.177E-03 

512 

34.7 

46 

40.9 

712 

0.125E-03 

1024 

26.6 

87 

32.5 

1042 

0.1.34E-03 

2048 

18.7 

158 

22.5 

1065 

0.597E-03 

4096 

12.2 

281 

14.2 

1127 

0.529E-03 

8192 

7.2 

502 

8.4 

1104 

0.461E-03 

64 

47.6 

9 

58.2 

123 

0.230E-02 

128 

60.7 

25 

77.3 

479 

0.180E-03 

256 

64.1 

59 

81.2 

1204 

0.124E-03 

512 

62.5 

128 

76.3 

2492 

0.125E-04 

1024 

58.8 

267 

69.3 

4672 

0.862E-05 

Chapter  5 

Generalizations  and 
Applications 


A  new  class  of  bases  for  £^[0, 1]  has  been  constructed  in  which  a  variety  of  integral 
operators  are  represented  as  sparse  matrices.  The  inverses  of  these  matrices  are 
also  sparse,  a  fact  which  enables  the  corresponding  integral  equations  to  be  sob'ed 
rapidly. 

Vector-j.pace  analogues  of  the  bases  were  also  developed,  and  the  latter  bases 
appear  preferable,  due  to  their  flexibility,  for  the  solution  of  a  variety  of  problems. 

To  pose  integral  equations  as  finite-dimensional  numerical  problems,  the  Ny- 
strom,  or  quadrature,  method  was  used.  Quadrature  rules  based  on  the  trape¬ 
zoidal  rule,  corrected  to  yield  high-order  convergence,  were  devc'  vped  for  several 
types  of  singularities.  The  combination  of  these  quadratures  and  the  vector- 
space  wavelet  bases  yields  a  fast  algorithm  for  the  solution  of  second-kind  in¬ 
tegral  equations.  We  have  proven  a  time  complexity  of  order  O(nlog^n),  but 
observed  order  0(n)  performance  in  practice,  where  n  is  the  number  of  points  in 
the  discretization.  This  cost  should  be  contr2isted  with  a  cost  of  order  O(n^)  for 
direct  application  of  a  dense  matrix,  and  order  O(n^)  for  direct  inversion. 

A  number  of  limitations  exist  in  the  procedures  described  above.  These  re¬ 
strictions  may  be  categorized  eis  “software  limitations”  and  “research  questions”. 
We  discuss  software  limitations  first. 


Software  Limitations 

In  both  the  function-space  and  vector-space  settings,  we  have  assumed  that  the 
size  of  the  problem  n  has  the  form  n  =  2^k  for  some  1.  This  restriction  is  not 
fundamental;  it  merely  simplifies  the  software.  For  the  multi-wavelet  ba.ses,  one 
could  expand  “deeper”  in  some  parts  of  the  interval  of  definition  than  others,  to 
efficiently  resolve  the  the  right-hand-side.  A  similar  effect  can  be  obtained  in  the 
discrete  Ccise  from  non-uniform  spacing  of  the  discreti:;ation  points. 
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A  second  software  restriction  is  the  assumption  of  only  diagonal  singularities. 
This  case  is  an  important  one  in  practice,  but  in  certain  situations  one  may 
encounter  singularities  or  near-singularities  off  the  main  diagonal.  The  scheme 
described  in  Section  4.3.2  for  transformation  of  a  matrix  to  wavelet  bases  can 
be  readily  revised  to  an  adaptive  scheme,  which  works  as  follows:  an  m  x  m 
submatrix  A  is  transformed  to  wavelet  coordinates,  under  the  assumption  that 
it  can  be  approximated  to  high  precision  in  each  direction  by  a  polynomial  of 
degree  less  than  k.  This  assumption  is  then  checked  by  dividing  A  into  four 
submatrices,  each  of  dimension  m/2  x  m/2,  transforming  each  submatrix,  and 
“glueing”  the  pieces  together.  If  the  results  from  the  two  computations  match 
(to  high  precision),  no  further  refinement  of  the  original  submatrix  is  needed. 
Otherwise,  the  procedure  is  repeated  recursively  on  the  m/2  x  m/2  submatrices. 
The  cost  of  this  adaptive  procedure  is  roughly  5  times  as  great  as  the  cost  of  a 
static  procedure  in  which  the  structure  of  the  singularities  is  known  a  priori. 


Research  Questions 


The  list  of  research  issues  is  of  course  much  longer.  One  of  the  most  pressing 
issues  is  the  generalization  to  two  and  three  dimensions.  Although  the  general¬ 
ization  of  the  multi-wavelet  bases  to  several  dimensions  is  described  in  Chapter  2, 
and  an  analogous  construction  creates  vector-space  bases  for  several  dimensions, 
quadratures  corresponding  to  those  of  Chapter  3  have  not  been  developed  for 
two  and  three-dimensional  problems. 

There  exist  integral  equations,  as  observed  in  the  numerical  examples  of  the 
previous  chapter,  in  which  the  inverse  operator’s  wavelet  representation  is  less 
sparse  than  that  of  the  original  integral  operator.  In  such  Ccises  it  would  be 
advantageous  to  avoid  inversion  of  the  operator  and  instead  obtain  a  sparse  LU- 
factorization.  This  wjis  attempted  with  the  btises  developed  in  this  thesis,  but 
for  the  operators  we  examined  the  LU  factors  have  substantial  fill-in.  We  believe, 
however,  that  a  somewhat  different  clciss  of  bases  can  be  constructed  with  the 
property  that  the  LU  factors  preserve  the  sparsity  of  the  operator  matrix. 

Another  question  is  whether  similar  “custom-constructed”  bcises  can  be  used 
to  create  sparse  representations  of  integral  operators  with  oscillatory  kernels. 
Initial  efforts  in  this  direction  for  a  limited  cleiss  of  such  operators,  in  particular  for 
Fourier  transforms  with  non-equispaced  points  and  frequencies,  appear  promising 
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Applications 

In  this  thesis  the  primary  application  of  our  new  wavelet  bases  has  been  the 
solution  of  second-kind  integral  equations.  The  bases  are  very  effective  for  the 
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fast  solution  of  a  wide  claiss  of  such  problems.  In  addition,  many  other  classes  of 
problems  can  be  solved  efficiently  using  these  techniques.  We  list  a  few  of  these 
problem  types. 

1.  Elliptic  partial  differential  equations  rewritten  as  integral  equations  by 
the  Lippman-Schwinger  method,  in  which  the  Green’s  functions  are  non- 
oscillatory. 

2.  Evolution  of  homogeneous  parabolic  PDEs  with  constant  or  periodic  bound¬ 
ary  conditions,  by  explicit  time  steps.  This  method  consists  of  repeated 
squarings  of  the  operator  for  a  single  time  step,  leading  to  an  order  0(n  log  t) 
algorithm  for  evolving  an  n-point  discretization  for  t  time  steps. 

3.  Evolution  of  general  parabolic  PDEs  by  implicit  time  steps,  in  which  the 
elliptic  problem  on  each  time  step  is  solved  in  wavelet  coordinates. 

4.  Evolution  of  hyperbolic  PDEs  by  a  method  of  operator  squaring  analogous 
to  the  scheme  proposed  for  homogeneous  parabolic  PDEs  above. 

5.  Problems  of  potential  theory  and  pseudo-differential  operators. 

6.  Signal  compression,  including  signals  of  seismic,  visual,  and  vocal  origin. 
There  is  also  reason  to  expect  that  analysis  of  such  compressed  data  will 
be  simpler  than  analysis  of  data  resulting  from  less  efficient  compression 
schemes. 

In  this  thesis  we  strayed  from  the  original  mathematical  definition  of  wavelets 
to  construct  classes  of  bases  tailored  for  numerical  computation.  The  basis  func¬ 
tions’  (or  vectors’)  principal  properties  of  local  support  and  vanishing  moments 
lead  to  sparse  representations  of  functions  and  operators  that  are  smooth  except 
at  a  finite  number  of  singularities.  There  is  little  doubt  that  other  bases  can  be 
constructed  along  similar  lines  to  possess  various  properties.  One  current  chal¬ 
lenge  is  the  construction  of  bases  suitable  for  the  efficient  representation  of  a 
variety  of  oscillatory  operators. 
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